Coverage for pySDC/projects/DAE/problems/WSCC9BusSystem.py: 100%
214 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1import numpy as np
2from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae
3from pySDC.core.Errors import ParameterError
6def WSCC9Bus():
7 """
8 Returns the Ybus for the power system.
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 = {}
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 1.0,
29 1.0,
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 0.0,
40 0.0,
41 0.0,
42 0.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 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 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 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 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 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 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 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 0.0,
170 0.0,
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 0.0,
189 0.0,
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 0.0,
208 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 0.0,
227 0.0,
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 0.0,
246 0.0,
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 0.0,
265 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 0.0,
284 0.0,
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 0.0,
303 0.0,
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 0.0,
322 0.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,
350 0.0,
351 0.0,
352 0.0,
353 0.0,
354 0.0,
355 0.0,
356 0.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,
373 0.0,
374 0.0,
375 0.0,
376 0.0,
377 0.0,
378 0.0,
379 0.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,
396 0.0,
397 0.0,
398 0.0,
399 0.0,
400 0.0,
401 0.0,
402 0.0,
403 ],
404 ]
405 )
407 return ppc_res
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 3.307378962025306 - 39.30888872611897j,
421 -1.942191248714727 + 10.51068205186793j,
422 0 + 0j,
423 0 + 0j,
424 0 + 0j,
425 -1.36518771331058 + 11.60409556313993j,
426 ],
427 [
428 0 + 0j,
429 0 + 0j,
430 0 + 0j,
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 ],
438 [
439 0 + 0j,
440 0 + 0j,
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 ],
449 [
450 0 + 0j,
451 0 + 0j,
452 0 + 0j,
453 0 + 0j,
454 0 + 0j,
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 -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 -1.36518771331058 + 11.60409556313993j,
476 0 + 0j,
477 0 + 0j,
478 0 + 0j,
479 -1.187604379291148 + 5.975134533308591j,
480 2.552792092601728 - 17.33823009644852j,
481 ],
482 ],
483 dtype=complex,
484 )
486 return ybus
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 3.307378962025306 - 39.30888872611897j,
500 -1.36518771331058 + 11.60409556313993j,
501 -1.942191248714727 + 10.51068205186793j,
502 0 + 0j,
503 0 + 0j,
504 0 + 0j,
505 ],
506 [
507 0 + 0j,
508 0 + 0j,
509 0 + 0j,
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 ],
517 [
518 0 + 0j,
519 0 + 0j,
520 0 + 0j,
521 -1.942191248714727 + 10.51068205186793j,
522 0 + 0j,
523 3.224200387138842 - 15.84092701422946j,
524 0 + 0j,
525 0 + 0j,
526 -1.282009138424115 + 5.588244962361526j,
527 ],
528 [
529 0 + 0j,
530 0 + 0j,
531 0 + 0j,
532 0 + 0j,
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 -1.617122473246136 + 13.69797859690844j,
547 2.772209954136233 - 23.30324902327162j,
548 -1.155087480890097 + 9.784270426363173j,
549 ],
550 [
551 0 + 0j,
552 0 + 0j,
553 0 + 17.06484641638225j,
554 0 + 0j,
555 0 + 0j,
556 -1.282009138424115 + 5.588244962361526j,
557 0 + 0j,
558 -1.155087480890097 + 9.784270426363173j,
559 2.437096619314212 - 32.15386180510696j,
560 ],
561 ],
562 dtype=complex,
563 )
565 return ybus
568# def get_YBus(ppc):
570# ppci = ext2int(ppc)
571# Ybus, yf, yt = makeYbus(ppci['baseMVA'],ppci['bus'],ppci['branch'])
573# return ppc['Ybus']()
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]_.
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.
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.
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 """
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()
774 self.baseMVA = self.mpc['baseMVA']
775 self.ws = 2 * np.pi * 60
776 self.ws_vector = self.ws * np.ones(self.m)
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 )
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 )
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 )
825 self.bus = self.mpc['bus']
826 self.branch = self.mpc['branch']
827 self.gen = self.mpc['gen']
828 self.YBus = get_initial_Ybus()
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()
836 # excitation limiter vmax
837 # self.vmax = 2.1
838 self.psv_max = 1.0
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
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!")
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]
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]
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]
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]
869 self.IC = list(zip(self.IC1, self.IC2, self.IC3, self.IC4, self.IC5, 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
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
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]
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, :]
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, :]
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, :]
913 # Calculate MH
914 self.MH = 2 * self.H / self.ws
916 # Represent QG as complex numbers
917 self.QG = self.QG.astype(complex)
919 # Convert VG0 and THG0 to complex phasors
920 self.Vphasor = self.VG0 * np.exp(1j * self.THG0)
922 # Calculate Iphasor
923 self.Iphasor = np.conj(np.divide(self.PG[:m] + 1j * self.QG[:m], self.Vphasor))
925 # Calculate E0
926 self.E0 = self.Vphasor + (self.Rs + 1j * self.Xq) * self.Iphasor
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)))
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
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 )
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
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
959 self.alpha = 2
960 self.beta = 2
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)
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)
969 self.t_switch = None
970 self.nswitches = 0
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)`.
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.
985 Returns
986 -------
987 f : dtype_f
988 The right-hand side of f (contains two components).
989 """
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 )
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 )
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]
1017 # line outage disturbance:
1018 if t >= 0.05:
1019 self.YBus = self.YBus_line6_8_outage
1021 self.Yang = np.angle(self.YBus)
1022 self.Yabs = np.abs(self.YBus)
1024 COI = np.sum(w * self.MH) / np.sum(self.MH)
1026 # Voltage-dependent active loads PL2, and voltage-dependent reactive loads QL2
1027 PL2 = np.array(self.PL)
1028 QL2 = np.array(self.QL)
1030 V = V.T
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]
1040 sum1 = np.sum(Vectorized_mag1 * np.cos(Vectorized_angle1), axis=1)
1041 sum2 = np.sum(Vectorized_mag1 * np.sin(Vectorized_angle1), axis=1)
1043 VG = V[: self.m]
1044 THG = TH[: self.m]
1045 Angle_diff = Delta - THG
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]
1054 sum3 = np.sum(Vectorized_mag2 * np.cos(Vectorized_angle2), axis=1)
1055 sum4 = np.sum(Vectorized_mag2 * np.sin(Vectorized_angle2), axis=1)
1057 # Initialise f
1058 f = self.dtype_f(self.init)
1060 t_switch = np.inf if self.t_switch is None else self.t_switch
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)
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
1126 eqs.append((1.0 / self.TCH) * (-TM + PSV) - dTM) # (10)
1127 eqs.append(
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]
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
1151 def u_exact(self, t):
1152 r"""
1153 Returns the initial conditions at time :math:`t=0`.
1155 Parameters
1156 ----------
1157 t : float
1158 Time of the initial conditions.
1160 Returns
1161 -------
1162 me : dtype_u
1163 Initial conditions.
1164 """
1165 assert t == 0, 'ERROR: u_exact only valid for t=0'
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
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.
1192 The state function for this problem is given by
1194 .. math::
1195 h(P_{SV,1}(t)) = P_{SV,1}(t) - P_{SV,1, max}
1197 for :math:`P_{SV,1,max}=1.0`.
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.
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 """
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
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
1229 def count_switches(self):
1230 """
1231 Setter to update the number of switches if one is found.
1232 """
1233 self.nswitches += 1