From 9263699909a658af38008da4160f7efe74c6b05c Mon Sep 17 00:00:00 2001 From: lukelowry Date: Sun, 31 May 2026 16:49:37 -0500 Subject: [PATCH 1/3] Document ESDC1A exciter model --- .../PhasorDynamics/Exciter/ESDC1A/README.md | 216 ++++++++++++++++++ .../Model/PhasorDynamics/Exciter/README.md | 1 + .../Figures/PhasorDynamics/ESDC1A_diagram.png | Bin 0 -> 58305 bytes 3 files changed, 217 insertions(+) create mode 100644 GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md create mode 100644 docs/Figures/PhasorDynamics/ESDC1A_diagram.png diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md new file mode 100644 index 000000000..e70e9c74e --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md @@ -0,0 +1,216 @@ +# **Exciter Model (ESDC1A)** + +## Block Diagram + +Standard model of the ESDC1A Exciter. + +Notes: +- Internal voltage signals are on model base unless otherwise stated. +- The source diagram labels the optional multiplier input as `Speed`; GridKit uses machine speed deviation, so the enabled multiplier is $1+\omega$. +- The PowerWorld parameter table names the UEL selector `UEL`; the block diagram routes `UEL >= 2` through the input-error summing junction and `UEL < 2` through the high-value gate. + +
+ + + Figure 1: Exciter ESDC1A model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +
+ +## Model Parameters + +Symbol | Units | Description | Typical Value | Note +------------------------------------|----------|---------------------------------------------------------|---------------|------ +$T_R$ | [sec] | Transducer time constant | TBD | Block name: `Tr`; if zero, $V_C$ is algebraic +$K_A$ | [p.u.] | Voltage-regulator gain | TBD | Block name: `Ka` +$T_A$ | [sec] | Voltage-regulator time constant | TBD | Block name: `Ta` +$T_B$ | [sec] | Lag time constant for the voltage-regulator input lead-lag block | TBD | Block name: `Tb` +$T_C$ | [sec] | Lead time constant for the voltage-regulator input lead-lag block | TBD | Block name: `Tc` +$V_R^{\max}$ | [p.u.] | Maximum voltage-regulator output | TBD | Block name: `Vrmax` +$V_R^{\min}$ | [p.u.] | Minimum voltage-regulator output | TBD | Block name: `Vrmin` +$K_E$ | [p.u.] | Exciter field-resistance line-slope margin | TBD | Block name: `Ke` +$T_E$ | [sec] | Exciter time constant | TBD | Block name: `Te` +$K_F$ | [p.u.] | Stabilizing feedback gain | TBD | Block name: `Kf` +$T_{F1}$ | [sec] | Feedback lead time constant | TBD | Block name: `Tf1` +$s_{\mathrm{spd}}$ | [binary] | Speed multiplier flag | TBD | Block name: `Spdmlt`; nonzero source values are enabled +$E_1$ | [p.u.] | First saturation voltage point | TBD | Block name: `E1` +$S_E(E_1)$ | [p.u.] | Saturation value at $E_1$ | TBD | Block name: `Se1` +$E_2$ | [p.u.] | Second saturation voltage point | TBD | Block name: `E2` +$S_E(E_2)$ | [p.u.] | Saturation value at $E_2$ | TBD | Block name: `Se2` +$I_{\mathrm{UEL}}$ | [integer] | Under-excitation limiter input-location selector | TBD | Block name: `UEL`; 0/1 = HV gate input, 2/3 = input-error summing junction +$s_{\mathrm{lim}}$ | [binary] | Exciter feedback lower-limit flag | TBD | Block name: `exclim`; nonzero enables zero lower limit on $V_{\mathrm{FE}}$ + +### Parameter Validation + +Invalid ESDC1A parameter sets are rejected by the following checks. Source data may apply PowerWorld-style autocorrections before these equations are evaluated. + +The required checks are: + +```math +\begin{aligned} + &K_A > 0 \\ + &T_R \ge 0,\quad T_A > 0,\quad T_B > 0,\quad T_C \ge 0,\quad T_E > 0,\quad T_{F1} \ge 0 \\ + &V_R^{\min} \le V_R^{\max} \\ + &s_{\mathrm{spd}}, s_{\mathrm{lim}} \in \{0,1\} \\ + &I_{\mathrm{UEL}} \in \{0,1,2,3\} +\end{aligned} +``` + +The saturation points are either disabled together, + +```math +\begin{aligned} + S_E(E_1) = 0,\quad S_E(E_2) = 0 +\end{aligned} +``` + +or define a valid two-point quadratic saturation fit: + +```math +\begin{aligned} + &E_1 > 0,\quad E_2 > 0,\quad E_1 \ne E_2 \\ + &S_E(E_1) > 0,\quad S_E(E_2) > 0,\quad S_E(E_1) \ne S_E(E_2) +\end{aligned} +``` + +### Model Derived Parameters + +The UEL routing flag and off-mode flag complements are: + +```math +\begin{aligned} + s_{\mathrm{UEL}} &= + \begin{cases} + 1 & I_{\mathrm{UEL}} \ge 2 \\ + 0 & I_{\mathrm{UEL}} < 2 + \end{cases} \\ + s_{\mathrm{UEL}}^{\mathrm{off}} &= 1 - s_{\mathrm{UEL}} \\ + s_{\mathrm{lim}}^{\mathrm{off}} &= 1 - s_{\mathrm{lim}} +\end{aligned} +``` + +The saturation curve is fitted from the two supplied saturation points. If both saturation factors are zero, use $S_A=0$ and $S_B=0$. Otherwise: + +```math +\begin{aligned} + C &= \sqrt{\dfrac{S_E(E_2)}{S_E(E_1)}} \\ + S_A &= \dfrac{C E_1 - E_2}{C - 1} \\ + S_B &= \dfrac{S_E(E_1)}{(E_1 - S_A)^2} +\end{aligned} +``` + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------------------|--------|---------------------------------------------------------|------ +$E_{\mathrm{fd}}'$ | [p.u.] | Field-voltage state before optional speed multiplier | State 1 in Fig. 1; source label: `EFD` +$V_C$ | [p.u.] | Sensed compensated voltage | State 2 in Fig. 1; source label: `Sensed Vt` +$V_R$ | [p.u.] | Voltage-regulator output | State 3 in Fig. 1; source label: `VR` +$V_F$ | [p.u.] | Stabilizing feedback washout output | State 4 in Fig. 1; source label: `VF`; algebraic when $T_{F1}=0$ +$x_{\mathrm{LL}}$ | [p.u.] | Lead-lag block state | State 5 in Fig. 1; source label: `Lead-Lag` + +#### Algebraic + +Symbol | Units | Description | Note +------------------------------------|----------|---------------------------------------------------------|------ +$e_V$ | [p.u.] | Voltage-regulator input error before lead-lag block | Includes voltage reference, sensed voltage, stabilizing feedback, stabilizer input, and selected UEL input +$V_{\mathrm{LL}}$ | [p.u.] | Lead-lag block output | Input to high-value gate +$V_{\mathrm{HV}}$ | [p.u.] | High-value gate output | Input to voltage regulator +$S_E$ | [p.u.] | Saturation coefficient evaluated at $E_{\mathrm{fd}}'$ | Uses derived saturation curve +$V_{\mathrm{FE}}$ | [p.u.] | Exciter feedback signal after optional lower limit | Lower limited at zero when $s_{\mathrm{lim}}=1$ +$E_{\mathrm{fd}}$ | [p.u.] | Field-voltage output | Output after optional speed multiplier + +### External Variables + +#### Differential + +None. + +#### Algebraic + +Symbol | Units | Description | Note +------------------------------------|--------|---------------------------------------------------------|------ +$E_C$ | [p.u.] | Compensated terminal voltage magnitude | Source label: `EC` +$V_{\mathrm{ref}}$ | [p.u.] | Voltage-control reference | Source label: `VREF` +$V_S$ | [p.u.] | Stabilizer input signal | Source label: `VS`; optional, defaults to zero +$V_{\mathrm{UEL}}$ | [p.u.] | Under-excitation limiter input | Source label: `VUEL`; optional, defaults to zero +$\omega$ | [p.u.] | Machine speed deviation | Source label: `Speed`; optional when $s_{\mathrm{spd}}=0$ + +## Model Equations + +### Differential Equations + +```math +\begin{aligned} + 0 &= -T_R\dot V_C - V_C + E_C \\ + 0 &= -T_B\dot x_{\mathrm{LL}} - x_{\mathrm{LL}} + e_V \\ + 0 &= -T_A\dot V_R + + \text{antiwindup}\!\left( + V_R, + -V_R + K_A V_{\mathrm{HV}}, + V_R^{\min}, + V_R^{\max} + \right) \\ + 0 &= -T_E\dot E_{\mathrm{fd}}' + V_R - V_{\mathrm{FE}} \\ + 0 &= -T_E T_{F1}\dot V_F - T_E V_F + K_F(V_R - V_{\mathrm{FE}}) +\end{aligned} +``` + +CommonMath defines the [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) target and smooth approximation. + +### Algebraic Equations + +The algebraic targets use CommonMath helper notation where applicable: + +```math +\begin{aligned} + 0 &= -e_V + V_{\mathrm{ref}} + V_S + s_{\mathrm{UEL}}V_{\mathrm{UEL}} - V_C - V_F \\ + 0 &= -T_B V_{\mathrm{LL}} + T_C e_V + (T_B - T_C)x_{\mathrm{LL}} \\ + 0 &= -V_{\mathrm{HV}} + + s_{\mathrm{UEL}}V_{\mathrm{LL}} + + s_{\mathrm{UEL}}^{\mathrm{off}}\text{max}(V_{\mathrm{LL}}, V_{\mathrm{UEL}}) \\ + 0 &= -S_E + S_B\,q(E_{\mathrm{fd}}' - S_A) \\ + 0 &= -V_{\mathrm{FE}} + + s_{\mathrm{lim}}^{\mathrm{off}}(K_E + S_E)E_{\mathrm{fd}}' + + s_{\mathrm{lim}}\rho\!\left((K_E + S_E)E_{\mathrm{fd}}'\right) \\ + 0 &= -E_{\mathrm{fd}} + \left(1 + s_{\mathrm{spd}}\omega\right)E_{\mathrm{fd}}' +\end{aligned} +``` + +CommonMath defines the helper targets and smooth approximations for [max](../../../../CommonMath.md#derived-functions) and the primitives [ramp and quadratic ramp](../../../../CommonMath.md#primitives) $\rho$ and $q$. + +## Initialization + +The machine initializes $E_{\mathrm{fd}}$ first. For a standard unsaturated start, ESDC1A reads that value along with attached $\omega$, $E_C$, $V_S$, and $V_{\mathrm{UEL}}$, sets all internal derivatives to zero, and evaluates: + +```math +\begin{aligned} + E_{\mathrm{fd},0}' &= \dfrac{E_{\mathrm{fd},0}}{1 + s_{\mathrm{spd}}\omega_0} \\ + S_{E,0} &= S_B\,q(E_{\mathrm{fd},0}' - S_A) \\ + V_{\mathrm{FE},0} + &= s_{\mathrm{lim}}^{\mathrm{off}}(K_E + S_{E,0})E_{\mathrm{fd},0}' + + s_{\mathrm{lim}}\rho\!\left((K_E + S_{E,0})E_{\mathrm{fd},0}'\right) \\ + V_{R,0} &= V_{\mathrm{FE},0} \\ + V_{\mathrm{HV},0} &= \dfrac{V_{R,0}}{K_A} \\ + V_{C,0} &= E_{C,0} \\ + V_{F,0} &= 0 \\ + x_{\mathrm{LL},0} &= V_{\mathrm{LL},0} = e_{V,0} = V_{\mathrm{HV},0} \\ + V_{\mathrm{ref},0} + &= e_{V,0} + V_{C,0} + V_{F,0} - V_{S,0} - s_{\mathrm{UEL}}V_{\mathrm{UEL},0} +\end{aligned} +``` + +This closed-form start requires $1 + s_{\mathrm{spd}}\omega_0 \ne 0$, $V_R^{\min} \le V_{R,0} \le V_R^{\max}$, and, when $s_{\mathrm{UEL}}=0$, $V_{\mathrm{HV},0} \ge V_{\mathrm{UEL},0}$. Saturated voltage-regulator starts and active high-value-gate starts are outside these closed-form equations. + +## Model Outputs + +Output | Units | Description | Note +----------------|--------|-------------------------------------|------ +`efd` | [p.u.] | Field-voltage output | $E_{\mathrm{fd}}$ +`vc` | [p.u.] | Sensed compensated voltage | $V_C$ +`vr` | [p.u.] | Voltage-regulator output | $V_R$ +`vf` | [p.u.] | Stabilizing feedback state | $V_F$ +`se` | [p.u.] | Saturation coefficient | $S_E$ +`vfe` | [p.u.] | Exciter feedback signal | $V_{\mathrm{FE}}$ diff --git a/GridKit/Model/PhasorDynamics/Exciter/README.md b/GridKit/Model/PhasorDynamics/Exciter/README.md index 81144066b..f204ca63f 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/README.md @@ -14,4 +14,5 @@ device internal voltage. There are a few standard Exciter models - IEEE Type 1 Excitation Model (See [IEEET1](IEEET1/README.md)) - IEEE DC1 Excitation Model (See [EXDC1](EXDC1/README.md)) +- IEEE DC1A Excitation Model (See [ESDC1A](ESDC1A/README.md)) - Simplified Excitation System Model (See [SEXS-PTI](SEXS-PTI/README.md)) diff --git a/docs/Figures/PhasorDynamics/ESDC1A_diagram.png b/docs/Figures/PhasorDynamics/ESDC1A_diagram.png new file mode 100644 index 0000000000000000000000000000000000000000..c77a100a65fa9450a6c12dbecf67e21b8ce45278 GIT binary patch literal 58305 zcmZ_0by$?$7dAS8fC7VbmxIzJ(kT-9l7ggkcXyYxlnMw)h!O&l(p}OZ(j7xcmvnvm z!QXezb^bVWi8{l~^XypbUTfX=+7t0o@i`oe0t*I%!T)<9r2>Q9p@hLuYwluzSFU9k zz!MmXv&!>lu%ZE~b@1elxrBlQ3|1P8eQAgep5JqLq2&yN5x7750T zD{V;|^jH4W;`Wxs$!5Q%vB*2}eF#lwgO^H1VuUWy6s1wIJR^%KklL)%aZ0T4*#IIh zG>Zv6O1|E1f9h&F(yjC%K}67W#07TkyzHHZ`Tmi68INvfpKDMOm?t4W{3btN@7HI7 z+M-?sc8&f+_$wAVy5jdYBt+88cLh$T3zWByqksQ@K7?X3eOxpU^?-1t9qBkx8EwO- zJ+er8M7?}b%GO)w??`!r)Jj~jr6)XD*JljnnW&^b)sz_gtYtH8_d4@wy%<=q6tkg= zZU6gV9X&GIy|KG)QSe zjeQwouyf9Q}>?S4#>a=|J1*7Yy&aR-CwS&|;9yWo%3+M@*Nth8fwWT6%bh z1Xy$7M} zrx72{tRhh(3)Y=-Rf!%Hh9F72KZGg2!4E_^6Z}9ypkm4OnIHkf!iqXaC0cgwFe2h! zO*j39=(3ZbMat@E;4M;YVB+xMjkfUgH?tlf9A4&2vI(un2tTG*v7$-JTsYkHL-+1yA>qG`gGYSbdCge3aXzK8PnpprH0ldBKq_kt?N zT}AhNi^DDjVvs~cel9-4#yu#RnS;M?{2splZ@1wQG?#z3Zg$L7dNUGuQcQXi*&H~* zk=pOINi5CG)XpT79coBN5$*LvvN_3ip{e0i!^Q6^ZPtcwlUQWLm3rTUYSL!RTfUP& zICNdH_+FL6m+I&|emYE?OCoo&2ve=OtB*S_sab!r7W`v!EL+RicHOSm`L1@8EK~H@ zarQs83T-c4FL*0m?h5;yZ0j?IpLB@RWBRnv?~6t9^_JCN#EHu}ir?ZQzlgxI+imPd zocV}x@)IaqUMq@zoutFV3z%>n*WK?FzaF~KcF~g#)7S}*94NVq(G|;Rx9d;k|6J5U z_>6B!u-Ct>A=X7GdU4TkSMUS9UQ|~Q-V>qI@rw-YxXf`psd^HvN4s@C^ZEH+fiaqF zp%D=gS!=lSn}oYzrdsG3sh3)jhD4J#B&C{#e}w-HrWe>$RwH8&{F|z*p7gD=@-A>; zAvV07D$yW6!5|v81;k5q^W=#!i|BNx|Cx=r%rzf`~Jc9=G@4-v-gfhIxPfz!Y3Q_`0#zI+j*aco-hSrHg7^x=d~L!#UJ0Z)7Rzt=#-lYBh! z^nWkIJ5u0MrX^;QSIa5z;aBfZuhCA7+R#o&Y^prgrB5e*M{#ipsQ$u*ZH~{a+WP-l zz$3}8Hec=mw>PoyNHYCT}JYd##)p^iz?{#f`i5p+KO{(~S^H-vmng5+FFt=0WsBU} zWT8jpExWI`scgt;+nN{fcKlq*>uZ}UrUW$|E}xfPMp7=TvA6ZzZba6>$H}Lv*=O`o zC0FB|hDSZh_7X2LQS5re1A=U)@h$$SMi7sq?fi;lM;aek#cIwEJpE@7YM}tkSt1_5 z?&;@NsqPsVICNo4`3!k6@h}Y)ug0>qUdBtq0{+&H9#{>YyB3_y-ENKRJ+SFke@nAT zJ#XgA)+2&k%it&TuYP^Yd4{Mn>^;@AaBB;_*llrg7rwQJ4I2nlEDRj(ju*d1gGraO zbd42jZ@cDw)2!JFJ+iiMTRP<6NS3%YW2v-YFeagJ+=IU(566ld6!hTL8q4>osf^Sm z$+OESe_PudtN*hAQpf^u6HNGRJoZxitn!{}-rm>nf4@$Geyw=Dh4Fa$O^K-w&U7wA z^61OYnc$Y2>c;yuq_Wxs(0r5`Z^Os+R%q7W8`py5YkhRVy1)?L2*|K&Up>qiTw9Sr zoWF-mt%kN63@b6M1_moRt&eiN%IBTn77MkmaNe81ZhFsepNdTIE#1|t>*-fGxh$S< zKc2slu)}>g$G;dgCMn$#y8i7yn~aopaFK<^HHppk?Pa<+U}|q9F~hKSuzz25o9&ch zD+>6C*CgNil@Hk25tzgzoRa1;XIv)X(46~?DR1~Wyg(zy=%QGkj&+lXXrZ>Pnz;0HSGiBMYzVRa=aY#;2!O>+5 z?<2fnBa-O({x9Yl0c+;ftWmUUVCj$igBg{i_XaHO3%j3S=~ai~-ho-YPmL_puX=5F zYyMf>02nt=!)W*hQ!x&VpdmQ@xLO#rKV!Z zn{)sRn=@nzUV= z7y!qgx%#u3B_9qMr1a*4-Q8hvpQi@YcOO@*ChkAhDM%yq&<8++RkTN2ycOkpGF4%a zs_t<;Zf|t4dqc7&(BjY}(8 zy!OVbJR#JA8=~W2VuVF*kZ_fM_5*MV-IOI3004Iixs3%c2IU4QCT^Y76Qul7?2qYx zQ)1t4aYAPwco(?OKRl%d1HFA6f_YHayv?xnmb2S<9bA5F@DHx>d8xr>|9gzQgs%4Q z0q+18!V(cbBB^Y;PM?v7nGMdsvBQqzQzaWN>)iL2j&BlEyIrTk%zFO?G2O7_|Mxoo z|IvlmgX@sT@9%Of_7oI1PHup8eN0bb)2srsuWt$lOYtVK<#X0H21XbczSZQs@RrD5QXcwt)c?MSlW5ckLfm$1dtMdKC*!p1k%vIQ_~IrG0~^(!IdFn zFVe3tTboB3SE~zS9zG`_ZrJUoOCJZgrI$&+#j6%{n`dlDfjlBAL#c_K|D^u@{ekvo zzq;OjRTPk$zsFrlW~YS6nAPWKdPHwoUf3q=)Mv*&`8J58A5 zR%-^(?r7-f^Vqa0`Ef3qsfWuGLtd5T*8$7`8$fb;8ZY(OT&SrzN1O(KJJ|O`eMhwG z4zIh~>DS5~vE|SWKn7D&W9)kq=-qJQbTu{2b+wtatYB3JR#3CSi3kRP3#Z)Y*Z(q0 zW&JM zgz$=2{zdfJzr<+yJMxwEN8Y1mxD9GkE+ z%5Bs+{qWRP?f2Jch9{fa*X?s`-J`Z9k!rQbqG<85b?xbCtp$< zyo`V$2X{0Y{YjJkbhAJ6-4npN)#IW6ie!EbyvFHDvY=v~|EBPzM>4^c^Kf*Xl;R!L{>S zKL<;GP8-}6cd%Epyh@0HW3wK6_`P~!9t7B{EH8zx6nGd9|qQg2P*#*e2k|FyRr}Dy_Ry%Yfr+ z)n(uF=Z)UX54o{>`t#`K-9-WDgvewh0Tjp?a1?F9sFrxspPCM9$9o?zUIzk#Y-7?A z-%|O#SJqD-$$wjWgqh9){1w_Pgx~~>c_AP!GkC(TZh9yqJw*y&F65A04ih3abARtZ z0b$k+1cRc{f&C!Zl88M>2Iw3+;B$g$)dGwhA`t9Ma1U^!<#&j_nn^9kl^mJ@`K^ySENzd zfN+b0b4R}J=tsMkKXp-w0on{s2(Su}B$;uB;}KEq5mwX&?F<%=Ah{=giQw?P>_2Co z_{D59EdwY?Loz_|99Z3gXHHx&9C)d@dBgVr#m1(@;HF^x{9(=@cx=O-@b0lh8_fc? zOZXzeN+lsmNag~1FucF*8)rmqc+n=PtOIYWG^;w77lwBJcuT(mw!-LuyGjQJfVwG| z9c4JeA0REtysy2Z`G>fO?e-5Fh9SSd34*l?ojg8ZJw4w(G5Xt@ik+zD?_XEX3|pAY zO=q35o9->{N{#R3Q4okhQy5-WyG>}=t9-2pHhn`Lgjl+pRhFb+YtJNi-BDhXkk{te z5N+-cz0J;n&7W6PL8W4!doU(|K4J4phcLU-syeW_wAxI$ z!jxb3w}vYFTpfEmu5Amuld;~MS$(?d9$rxpnbz||&(LuZBuV$+8S^hD5AU{Xv_x?uy713nQ4 zd*#R7>j7FE3yPc^N`m_9Zr85~U?reBpBqU9xq=}m;7Ck09LS`%ArlKhr>tza#8SVy z+V(j9w$_7zFXh%oIB=t+H%c%|V8f+A2-u_d-A~oj{_jw_TQK}{bxpXRQapQbIBpZL zE7n>W_QCh%uA5O_2S}DV^Pn^}5*5b}hi@bX@)XSylt#MNPc>l+G4^vog3kBj`B7kn z=NHqD&38UmZ*8P8D}En?w_wDNX6CbSC%>&vfwJMe^<@J|c3I-^Emi13yp(y&&^m#i zAKIOUYjKfuu`?%R11yjGg4{?#?gXf^JSv844g`i5z#@lrcFb5sE=HUKLui})E*bZ_4EbunVt5J^W*NMkBy`g+P(51z#y;dI-5IzUz}Qu#IXekJ@1J}#z zWFNhm8r|oOJyyZT%y%Aqd~)*A{tDdx=j3mm)HE_11I%jIQ@z__l@}9ugIxF&d=Fj_ zmj-;GbbCbu`VcUJD2K^W(#6fH`C!ObtD^7H$v4%{PMfVcmAX%dF#6pY6g#5gs5@T@ zI-2FOX@h*XEh|n+70}yYWu?9-+edUFOW5S*h#aa;>{&z2ptYf#uEdX!dKT+BQUhj} zSaZ1%g%yOGpy{zR1Wj{!72&)}vt-JSQ$2}OW<^s)tx5DQywxHlRFwr?TDiV8Ct03> zhwGaq#k=l@RGktm2}W0>VD+p7F(~s)A!}5>J{)y&&I2oi1FTW=k1cSf$g0zF+Dj>k zp7o@v-U6YRj~nmOP!U)NOoZ(;hbCvQS@}yf;{5!c0?5i!fQMLplMUsD3{Latq{?`7 zem;hG{?bP|4IKnKa1G=0mF>^T_vf9hpT7CN`w&qONR5VI>1c|@L+RnFw0$BM39d5* z7^w{b1lrO;0q4*BrY35DEz)!C-}_6PU;LcR5E12~b4;9W1h#=o{d$O5VDu5bn(Q4s zn`zvByArmj7CC5ItYBKkx~nAbF6Vq>(J`nldq&;z9)8Kv1$vL|v|a-SB1=UIuf_4x zudPtP!Ac<>@I(7rffEcLQO?ny+Zi>41qe5KV$;;v8%s1BrfU03;KV~Ky>~Rd+d+RK zx>^x+vEjmo^`BG)Q^G7=zUz|$K&*Im1GQyFiW;oN1{v7rPr^P+2+o(cmVr#IGFXhD z-QRc;{u$p7Tj}ye%ikpH#xu|6@7uq++ts>DA`hgbBxFtbxu8ph3rYMVp0Z@OJl53e zORkt}Y6`X*d(;2rWO4#A!4^DN(#mnYE6bxJs6!S+u+!^T8fCVWfC-Pxg0N1RXrdCX zJARhNfXGUJ5aS+d$VFJ)3#onzkbP_{4%z`+bSS6rK*)g5kp|L%EIvGBA+G|m+6s}M ztW(!g#)jJqxmirJA(846+%I{$JQN{FfrY{127$LhwP}z}H_J*u{=sNN02c^jSl00S z_a&8y?BlD*X=rfan&CK+QV3-p@0gEv!}A>Cd|zSz@ujj`&?6A<9Lc0fdkIJX#NHf+ z<#4PYsrwZ-ob;J*(Gq>PUjA(rYPgKD3~!6mP>y1be5XMqZO5>|@=cZHKHYt~3~17T zde0a%;r4``$uqnkyf~IzvX{7I2-@l7dRw2x^!cJlWwRZbt-@T)p47Xdrp7)I6CVD8 zA;!?KlWao=!pLdUpG69~p{sCY)4U(1rM2%+tqjNa_HP^!c|%g$7z)FYW$-&g9M=jA z?_UySTQFchE|jwivzP4$q`fDx%`}i2Zp)kydVx}tB7&gyYxcm2mw5z%sou>JIb*S!jS@VwV6*t z^nr|2RtG+rR|OI-YIc$+ln2ZWv!Sna36$%AklPHW`LMBPow5p#VB(*V7enfdyRae0D{h8Bn`-} z@-pyZJ0+Fcr4%qrcd@4JBMQNS)pu~cZHARvri6DG0ddTfHUp>yon>rNUxu1}ojJJ` z2;*Pe0rm;h%bQ4S|04$~%zck~!`Y7>>S72;1H`#gJFz<{vg zFXeovWOsQ(Dv7w9)I9K;;OWo0dX0z=u$`XcNp?_jF)sbcKZDrXj|k5L8tgN^RNLOh zF~@S`MuQb+dL43%lBI)30mz#^VzD|LVX7=4sw|@=6b_MDJUnt*YNJ`myq7x_ZfIbh z$cP{FGE(=aB}8^hJU1i4Q@S0gz7AXF*ZKu@(*9DQ#*w}b3;-kO=V0Lpjq%lYj zSIV+Uwx2(LVwolGmD^5Z6GUzQUU!jjdxI=ET<})}DuiZAj_2C%?@Ome*OZ})rlhCA zzY}ds%uv2e^I;ebIp1IOuRQ*DvbeXhd$4T zu)xg*MU#5F3%>Pd4Kcy?${zsdOn)@+9Nu=fA##psskAjO|Mk#WJ((UJ!&Ze9BEGrP zVTLVA#`KXC0^cHKEbfAER|`mtsE-m&feE;HYosN^M~G`(folFQA82CIAZ@AKjtu4b zP?u5itSxJhZX8{OD4VWL@(yjggV-<+b{Uh2k+$wQT1FF))r_Yqw@sbHZK*q%t+44) z;STd*b~H|1257G6k6EQfkqxi=D{TP~DGWc$eGB{^;4@<L@ih>$rCPjt%M5DK6Ktnt zgwkXn?HQ+%(g*Js z2m?hMLJ^O8$i3C~bPj%Rw;zMJdW%L_mG5cwMb=bV)J7P@m1-1bGHxmXfDNHtg@>B+ z$u}G;%8_h{&+>eDCiOzP{R42)@H|4`lVqXFZulCJz?I_e$OXGo7 zyRe5Lep7Qgngx{V2)c=IW$#)_p$uWX}7(gYbEc4sNpi z&jRxJ-j!!)@I3`fX2Rs5gN~Yi^vjKhAZ)-+HwB zB6xK>wX8eUNL(t1W)#h8htE}t%0vLm(6LEAG|)aqUukc@ZN%j=US%o(vzQ?1kANsY^><*n1C9N#&9_ItsPQ9r}644GvfP z(Wkd3O3{FZ4Iq_}kWkmVXlevC9@NDL8wgO4I=#|j!<(F!YWY|~^xZ8&3Zb__=88!Y z4oKH5JqG=n=gE^q6#tpr0r5*iXU2D>Cjs)ATf6O~a(-K1%QjuvFwPi+3qvh`AQk`o zal4U`5yCf>xKMV@@`Fubv9ZMQjwN5!v7kej6d<-hfzxCeOouevqJZ8UMrqkXr^6 z3S|&HS^mVO=74PA)$z<#uJ{#mQQeLZmFIH2W`h^6+oUO}=ljrj^ADLnt*!b)0b2)2 z#k)6vWPDX}8Flnhdh))MAL9HEt?;N~j5d+Vet*B5)QF{?_cpJN>!{jYD_XLa z?rj8paL#)-I{AtDQi6tFJZ)%r`1vO&rpAG2B_~3|SzE?=@heR14dkXbXBjtIT*BzC zxjANhcY6nB1o>70$4~mBEot^8+>9?yMBG+5;_S_nWOI!q;;n#Emp78I5n zj2UNr=6hq>_KRS_-)R7M;e8|UTq#>c+Lc1ew_G;#)#Wmycnj#$2>TK%Mp2+qoTxD_ zJ@O7OG%W0MzSZQK)Ts3aSESQevGo>d*!gorJWR-dDFl>1fXCS|;cPrb(0Fe>{^WPQ ze!v2T9vUC@==@tgR;ObRZZ<3+)@PYsv(kpVaX;Fg$44N9?3|56VVKOq;$0#JB_zK9 zO#sq6a{W(Xp1jH(Qh!TFaKBH=5C7fZ#lv?ki>Bl^+HE*7gW4@#XN~+?k`02x93TJ! z*j9b9Q@8&<<5EgbK{qqI1w)za0MeBuS;F(P^I2@Bi#)s-$LZenX{R`!NL}FTFW?GI!15qvt|mFp$VffX&4fREkrIaoPdS&8 zk`jrB1;TZNN*lfu9RN>FJz}ft%QMY-kH_hL=awRS&1eT!8xFSEe6VI9;U}hsvoQSB zsHn9I?y`ERR*|;{2lRsHRw+n996N8DFBCv~N;Wr@EWGq?nZx2Uo0)2SVYe-i5Jtub?`@OGK8Sp#&TV@iUs<@}`h0Vu+L2yS{K{pkW-T2e z=(0jth+p8XZDyCz2K9X`E$@BZzlXhnV+-d4+&82c0IhAg!uD2;uA%^x1sOyjL5hN@ zPvL@_)U$#x70pC+fEQjs*-HZhrfGw^9X~n#B>%U?gin4Z=P(x=86Jrdb%TVVrj#bO zsT$4yG$YMj5D4%<+nXaw3d6-U1T7TP5Elz|PXNirW117F6*S(9?sI(qlsqlGxR3ZK zBjOp`R|;3`fFM&FO*EkhKMQoO`BWcS`;vzO9p;xlQeITOgjIh$WnI0N-iq@KOw#PJ z0%4-*gZJGYq-%f*d5Lt{zh&53xp2nhHn|o0Z6N(0?8DF5Si9RFK&01`GXvcL!(F^~ z5}T(PF4pBr1wii{DbZWDr(oClyQu#C;V z_VhI>d}KC-!X|GF2x;x-II6tyMWcP6i9T(-3$N|}8BJNbSyTrycPT0MBzcC0g=gH! zCM1>_8+it8K#hCIkDu@eYY22EdaAXR0Jg`KGJdrAoEzK+H!_QhzIxDAwL}7t20+6v zwgoGPEv51wNX=cDLTsl&jtqiXVF^G!=5(Kf|LjSCL1yXwdVA<&RTP5-><+B)!mqYf zLm*N?O&(!oW@u_pr4&+235`G>h-+#mjW^EW+Ta=+)DbkoM@?#jC;bS>EVrWJ$A%!L zw;xfIUMmXUZwN73O+{9xl1t8mNDAQQY$sWPSv?@_L&ZmHDInq!lT=H@Fw(?F^f?? zmcn?d`st^2IPp8Hk zL;cT^{9D|cE_A6J+wmcH*??8r(kJumwbAYr^75HG-0l42fCa2#Z!4Czvl29@q^8_~ zeWHkJ;dz_oJm*IkA(d%KfnDiNo^$slAQ*&8$iC+1m)O=7y(j^o2vV`U>F`1f3<^M6 zyr+nzr-?=OyI*6@U$;1)4h{Op+TSb)7S@|mbI>;6Tkm*HSnE5Fw|N{k-GJJUDy3-L z)o&X6{OZX4OH8*7C$3N@g8D#N6!~OIg=t^CPWG#wI z{U%f4vpbzIxd6DyY^$pOvQ^;I&)v)HMEtfuHUe$tUf_Kn7#bQHV6h8U_{rMK*d`{~ z=KRI)gU;m$SlII~bB}Xq8~sNp*?L=jPi7(Jao6vp&WW(5^LphTJobhwiq~EAt+5PA zcj}v3=l<`fb64)j@U}0S8I|Ww;uCpWrdHo$?(qS{(H)QPHlY|Zoc&=*C?g}|FzY(r zUEnoF+PH^af3}(m7hDXkGnRkIbGa9$)t4;%MlcL+xOgwuV6h8f7^hq|5W)y?`ssajL& zXgKPOiG*8hU%mgFcJw+qXWjMhqTH8*T?=AqKwt*#kc@t-oC^ZnM-1E6Q_BxtgZ64c4Sr-W7{<|XiL{>07NcYu{K)#@_wx_OwJzN^;s?W93v zRG>_YvhyulYmau5%7Nx#-zk~;?5WnSz$=cD=Oc=UiC8Fb%ZTh9sy8X4)nQ0E%aFi@=#Z_*4+dG=3zDw$!G&D5-{k$G!(|nsNvBc_d~iI zsXjU7N7{`21j|yyl;_*mh48t+99)F%adS2`9{<8Z^zZ_aHD{8-;PoCdi>YQinyEbDnVnYCAs`$J#ot%CN)8qnV4a$1Dzrm*HhGL!fBb*P$8|Tsv{#oWYSFqF0lm-d z+xx40(3@sbg7bf#`992wF(t!D{M1yqboGghn`HH1h|=*afq!w685?jX2OH@RuqvaR zT{_bu20(!McHD}npb1m{i?{DpWk0X?*K4JzE&Ma#C`r-11F#V?oVvgF9uPCVDo6P} zl29q^9C9ig{=zpYOz#eSKJL|+sv6vvCL0Jr#v_LfC-*%qhB8Lq$M~i3lbnR~+VBETbPf0|) z#xbhI!slVqAo^XbZN&%36=~&2RQi7P?$Hhz@m#OQ6N&lv>eY*(ln&FjRj#%r;A96V zyf8qmMSHy$wlXUEIju9+?$x$^NxXz@wc6BfJ(8qQy9zhsYRi0ORLXWZwj~tAblfjr zzMyYk>1aQ)Bs)DR&St?aO!ry7CdQ-$dlGa>JNGj|BgSnTM&relfSf)QYiM)?lPgw%n+w`&6eaPiAN47D zoZ(j(u8ICU3ZHzhxoZSfzhQfn0T_JzSTLfy`@c(&RVQk=i-N5fIJ&B}l>yn5fn+dq z;}KaCtTH*4 zF8J)tHSrX^C)=L=-V1c$i%Zi-_KahKF~K%`(<23nPJz#xG}%bJM&~PS!HJc!Cm4WH zb;kTF1FzY`6l^CZ@pk~CQo8qU>tCn1mRK-oq(IEc1xQ{-=B?0GycVIyf8|VgRLkyfnjGn@q`Ez!meqn<%8fg5Gs9j*M&#A? zQBwS@Ok!MK4fSxVcbMwc(DTe**wuSr&i}Xf%%m7}z`#x(mLwXqfn-tS0@>M{ey3I8 z5H%p<`fu(*ekeA!z=sVT%Fv+-WTrM?-i!AvMQHA`2eC{Ai=l0>_t{y9>|fxM!yDWx z@>7iapCFA-+by3X7k6i(0iITYF85D0ZKvA>O(-yVB_(8blHNQ12MDWgFJ=9*9am@&(%1Idyjw^g;Vudg(9 z8r38)bs;>D9DZt*f{_g$u6^YT_Hgj8s;fpI1fS&tb((g#%lAPIWAA@49dy*6$rMTI zaJN_8fKT%1_d3tQnDJeTC#0&yfnJp5loY*=RQ64ya?ZI?gzYLJG{FF!EY`N32TMM( z7Y5Bh4hw;pp8hwh{KDl_K|uM{;aFo|5?HX`tawX+uxgM(1m2wQBOjq<`;FREovkyL z8nmDm+phK{qffiAcCt@wX2^(AM1x~75YGq9^>mAs-~CjD8iVq#gK%h*pC{_z3uIGk zVXY-PM}a}QOSNUm9&rgu0ROztaFsMNf?5sWi(@Bs*#5rbpW2;KGKiOSs?nL(<^92i z9si; zR}(tJ>5+#t*8{TxW+u?V3pqP%Sb^Oon;EUc)0Ube3Ct_ziK*zzNmC$!sHp)y)P#Gz zzx{{Hx(Y0yXcK*{v|*(QHuL+l+!1L4*s{!?(jzF=qqz5N16!TboD!B?eM@s|hTe`aL!WPWMFUYSi2XClS`Rzd`(Qe(h0UiZOv5b>M9S%$Y z>Ql1WixKki9wIa;S$$scXnLWPYw7ykN?Xh#Y^%e=V(Vp;lu3B>13&SlhjMD8=~GH+ zBpRi*`5pq=k3a&Qh~pssL@e~|@TT>mZltkLtKG57zm?6;+$iIZ0PgSgAMPCkKq8xG z#y`UCmAzeznj0F=1oHvV>-nh(2P(o+)y4YKtnE6$gal4LIb%(&^|Ok9cVKjydDaE2 zo^BILQo#YfCoPsKtLQZCf`vJj#QhA%MhKEjeqDgF{QnK{3ea*=NxBLEUcEBv|41%`b8p|6=y!4TiCYSN|h zX$v3Zc;DSsu<>&%^jy%c0Ykzi$kx4QQr6|bUUrxV2hOVLYP{Zj3<7mf2>98Z^}k=) zAhJjtfpe(D?}st;XAzyb1Ct7P3E#?x^VtssZd5Gpt!%VmQhUM{^RiO2?8)+gZL)Vu zdo6OAR<|E99D|D!Lj!nRh*ZX!5=y)KgCaB#s9mCrp2pfEkI#@<<~nD)LR zR~NX3&vx(48~T{j0F6FrfsMQ`Sl2~>_L`e0lF|9BHj~UMrTRnr{sf@?Qu{$C-NziQ zp@GVGyScX0T)rU#%O^9A*66qCf1F}&l#t9)PvK1_K~zV|*#5SAgcK?`Zt05&q-aUp zEkm#*l;9ZSXM%N?sf(FQ2@FRk5+W&rTEps)E(U8U+0wgKoVv9?=<*r>tHmRlt$Bxn z!jtMiKW%diiC2-{rD%4iFf!TIP9^%4gdOAal;= zus=N+H-~87%f9h!HK;GJd_bs|HBR$hr(1d zH+q6xpO_D{^N#BT>$**5XEsp8q?WOu-icnWHvZCaT5Lw3aa3hDyP-?I=@Ro{NjKU< zB*uh$vShyMSKh8KYXq8YcAfrE3(k3CRHimd0*eJH^#g#&>ev6NpcXRcsI>QcL8W>F+rJRu&8nH4F zdSzP6yiO4k;2CrB#zqyCF2RK__u8U>imx=VagV2Yn++YE!o^*-^fifuj+Pe)!TaO> zSL5g>N);K+0z1UVGME(S(MXcjznq320*unZ8Q{X)A*MdNyqG1D>%W`{g@=cjsBM+& zJe3VkQss%ri+Fk^z7jab;(b{s?+|bmmnh67Mz9y zEe_~O>lBUWqs2=t_Hz-3MzjY%qeKDBHGJTFyt48qTW$~F-OpFzGcYx^QlB%YGY$~j zo>6&kl{z`ANwN6-Qj}FD(rmQ6S`IGup>y*YzSDD3?sw%Sa(#+~G->6M7^_ag%yUrR z{`i21Jt8ev=h1#wLHU{^O|k7B8Q<`_6sFHGVXR?CSuE*SqKgR2L}}U5R}9(J$sxak zVkCi}?9Bm~TW}=0;?=VBA#n%}7GLwM{v1%UPWNN|-xm8{Z32FG1Ai^UUo4SbJM>nH z0%@44T1OB%HXuS}MFghkb3UB$lm5OxspFBV)SUmg@FZbNh#LR* zX0#JOxUgRxJU@S_XK1MTebbj1A*{6*hbZIEToQwXin8`%76fqL?NBuwymQK>C7kYM z^a1Y$-z@>Yf@15`$a=P!c42?BGDf*d`D-%HYkUj5kF&pk?C*(G9SR8)RT$j7NXg7)Y zQzb+g+qfc8C%HXO-^hPrgHlP4{SNVLEkb2`K_HVmY*&mMb-#NguoO;@8kM(OHSY%{ zvCv8a;PyOR?M`(ODOF-~=UA$iOz8mUYOo7uvR)*M0=cB=Vs_U7n3%zIr7e2wPyxZc zn(V7up$WQ1oe`Z*T7QwTBHgGVR`q5s939`I(K#bBj3TNA)Lza89`B~`Bv$$FaxA$NhP{i~d$Fzt zULwP{FaVcqv8WtA3IG0Jx2~2qpbEhFg80p$c&S@I1*f(7`tNLk^*2Z?mB^$|zUl_D zjt>c1erL-$Uha4H7Eo|(34DFgcF{w;o}lzL7-TS)faPp%y6kOQs2JCKuR*9B8!m+q zfK~}ekqbGgcC8TV60@osxVRAU&p4nJ)vQu{0v-9M&6r|z^E@ArW_u|)>uHkwa*uq= zz?WnwXG)B?aN)GoF)?uvM56D6Ncg9Idi`;tU?ty;R|i&B zgtOmi)i#9bg}nMX@FtnrNW!~@pRORpn#YxDC1eYvYc z9=va4tEXFP4Hsa%FLY9zGiZ^*FjlM|yV9G41T?+HW=;6r^Kd`w>Q>bFAj_S?$xO=t zAE7k?;ql(>=0}dqtsp85c8ad)-QY1fa(P2Lltc z6ly^LGTG0q$#Nh@IrcE;E|m3h1cMKv!5xI?6#M{dU#lZPlum9#Aj(}bikmrJ~gR_+LXo0&4Lcd()H{-pVmNNURYcOwCWk)r_p^5&wsc!T?#aIoERZVMle!Qm ziQ5vvm>%AAmfF;kD|Ym)78Sm?UEb@<2oj^o6AkCidu^nbMV;j4p8~&qW4_!n7e_io z;{KMqz~5j8q6W79DJpCMQUyUuvx#mmesea61Z7L2-4NN*r$;ogpmQ?JX}#gjR-Hp z` z4j4u-5J>sE=K~`tlfO2uUexS5XS12{{YVpNV#*x~w7e0l%2+BeCybVf;^vXVCyigh z4zJZVDSNkOKg|bLT*e*&(}(WD~tBJ7kkRLJ~q`k7Q*=$of4`b^VU}_#XH7IPO1h zhwJmXitBuz=j(ip=Xi`9o9V5DC~Rfrp9hs6T9j$>dh8ae?glk9YKYI=J^O?Ic$0r) ze9foKp9tUd|<*R+UW>%=ZG=Z zyck4ck^HV{{JOcc)~o*dop_&&ea>MltOxU-+ZMAk0K&vuWD}_*{mLMaWG?E~IA&zh zH2VWdnpcB+_BTv{@|0jRzulGxN`2qmV@)WeWB0hY6w?VKRkBMuMM?^!T zCdClY<|ArDZez$XoN>H&&Qm)8St&Aq`l7#EG>&`Vp0S~`|7It$Ud2UBOnleNt0wQ) z3wXEGTI3DymDUu0$|$_ec=_@I^Mg-UHZ8OnXn!f|aHgx`uL`QW7I3oIEmW9NL)+Ra z1n_zJyIa{Zz=LIWZmv8&+K2Y%!xcd3>1t-5%Gd4v+@%HpFryl3IJvwBoA7S3MA{Vb zg{MDH^0R_U&#-;v`1kGztzRjO&m_{gLy z_warLq9`J8V!$xo$C8#u?-791bm0*YI2**Wu6=9wR>!J`z1QmM8V>ivGaG_NzBh2& z#~l}zsT@h>s2sijBC@xpH2NKWQ1Z?0+ZsF8&$^w7hx;p)ApMc=dd&HG{BZu(Oj>Qq zgFTDi{HcI=Gr#&PlOwt-UkQ(C3XhrVq3!ZOTN7S(h>+8j4_7|{mLn2beM_lv3eL@< z*RB~B{}P{igpDi#z;NIC=)1ZN4={LRLu_y?bogEP$__kvL?MUx{k%g24r2cWl#i>#4mrn7j42edFEYDp-`zDxo+-Z zZsQ<%Q%LoLd9xuSkGA^#)kuKlsZTsHd=pjoRBP4wkMMA@SsNjbVO4SFVa}^4t`RyM zmgDa*23{~QOylgw;}Ifzj#v-B$0`>4?z>lu#|%dYr&EQNhZ=n&TOQDTe@wl9O*~%? z=X{=cbR`cFvOgC?aEb+|@cJnp{Swp!k5=`hr})umV%O2GrfkihPk1u#b4xX}&e_S< zN6xuaUwed6&>r+*kL!$%j-EgDebVQb%3hY#t=frH$9i!`sMegcHbhrAQnF}#L$Ed@ zFo4VNnD@k(OQzkfwDSLYAtxzBX%3K8SzWqDR z8J+9Ma>8BaqS;Exk=fd}%uvmCh61ixQy?x<=iC_3RO?JFGiva^p~dr4@+F5-n{Lh4 zqZx~G07&4oxW?@4+;ObaE`A{)KkN~b6Aww>?{N>4>go>I+}TBim4_YW+$;yV6n=`% znNotXn{!&=uZN@QtJTjDco`HbH;xiQ)5{_sMsM!&pKh?pU9CsmV71lCf2u~pLrr~EQSM@O%4y~?s<(ixB{9VzJYrJCBvy>`2X0O zmxOq8haRt_;b?2ec3aH-C~-9^s+mfT=`<(&z)%l+7VdHv=T1tF$42Sow_8_I=5ZbUJR(XSZXIy(R6h*NOTPt{&ADe(s=suZWWeYEApNIn ze@d42m`R>a?&Ngbt zA)-X^A$jz%b@&YrzKz|_ZWkPl0*dWX>dhsLbWRIPcY|hBE%j7lPOZ`n=DxgvdZnL$ z3lHpY>mAh!wp7W$3(5?#3V9Ho4jnuRuywJ zLFfy33egbR;;w|UxYJ)vR+bOUyxNw&@piybPJ6KTDAC_N(WOVtu!EyAeL--yXO=tn z{Z}3kgZk-BjnekO!*{8c;Z6ROLF@0W=Sknbr`d}!B$q?G^PM`UsRC&UmWgwm4KM{X zwTiwBbSmA4(N!(3GMI_}(Ps88FT=9(7JD4c;W}kFYC)ANEXk$dnQV~#sKO4nyT(zJ z23xMRlT3?}v8JH#pZyU>8{cEDL+PEG%mTcZl1O;~2&W;%H)Bw6E<{{G5)TCKSw1G+A zDT8h&pSPJhNo?o&;R>VYBA5L1y4fX0Hz`qRR`H{Sebki^4f=eS=5-p1gnF?&y!U+T ze|?rDIt@_EjtimG2jk&DBoShKTU}~Mc3gZO(qR9}Sx@a|gh!&k#TXJo!=Y4m_XZUg zKS%y4tl!@(&!xG3xrJ+Achgm(>*-~Cyep{QW2prH(5KfDwyB&;m|dA&5V8w%2pr~J zK;_1J?DKQu_8SecA4&IW7Zdx~+Ar9oEBUR-P3kU@!HD!o?L1E%Wru>Ra~kFzj1?Dy zEO`JyEQ7jD!6m|rFT(Ik} zj#&@CW{w$)&ZR)jL!1IEL!@_!pJgW!33`7Ghn-2dZ)#^(V3jQ@B&48i9yq{&WCDn=c0yx3nX8nae^#n#~ecIu{P ze#PuBzBgJrEl8lwAn3KZk}J9{S1%CxeZg22I=SW`E=fEl?I*V#B5b8M3? zwi?PiqZdip-|?;DDbiLJ+aAy)8fmG>KKg!K7dcyS;GNjR8$P07&K~w)&4T>!5t7dG zSjlX5T$JQm|H3z?+M>Xm^NAxqV*59PNh+W_KobvA-&V6LMe3&HeWql0CY4$5HD=(d zFPUqhsn;5{)flI#)UE}PFo`@jeDQDwv?asxm1^$1)v@VVaX>EZ1pv-;$Dg%Gg-BEWrizu;1*Ask{ho3IelbIaZJzifqcP$hF z)6BgRok@ySK-yQx&$u1+LnXiZ(;-Nu(YmTFcH$8Y_S%7<+s zA}GB(1^}4c%Vsj1znvpe^F|8izY?|9fAVMEj249i_~(pFC#l5O2%c~#w^K8R0$Tqm z1Shn1w-39Iz3Fr#=O(>-J)Zngl+H*@9Z~YJZ23I5eBpEa!F)#Zm3C&$bq^WxU^W)Z z`I*Zjb|$>|M)>UP&9*WFflVRoo7n{MOSLVB&Y9i1zBCFOg@K)mDkZKGMw;V94NQ); z6vta?$67Ac5T(zhdM;|t(BJt;8*nG+HC32w3%`Fk(&MFvk-6+hX$5_)uVdQ+Cm1US^LmgBuF z;xq2+3wyXOu?4b1@@y>EC`T%-1=WKl1nKpWbGO^kFX%RZUD@43gfKA^ewnsxJHkLh ztU~ftpG7W7P!Liz4*vKO|u8s%^p z0|`F^3gxMP!Sms0;oA+C!UK6JJ?i8Q@;K%KyOpugyqh9gQvOcHLvGNpada=Rp_2iL z%)MqJegtxky9_bCW13U*8C?4TnP$H=mzazOk0m^_b9u`Ws#4nnH9c=y*M$btab(>A zGuxpuz)CFg7~+(cmKNN|Di`j`Pzv2Cg%tJPL89jiN##MuI9xwN9Rb~{LgXve{iCXd z8L0rR`MNRV!x^yWQ=g@PrUoLk-DEurMBeWXV#aP&gT!T^1grd|4k#1!zPv1Mb1e!% zMt!_f&svoC5nyc4$T55Kz^!crVWpC?UIIJ}OdD@f1O#iU)1e?HjrrT}1#*0PTsZ%`2=D^rC~ej`m+-5w-&VQwFt z#5p8wmDo76-+p%shncVXJLQ^~){0~MWH^!pJI}@IeoVbP-5f0Td-xu(&#+mkRQiSZ zi+{O0wF%8(KT|T`C6bK1Y4Y!WR5@-6N+NT6NRDnmvrBs=}0zu_wC1I zb~Bc-{Yu~5N<@bf%j8U`J_+wdX~7MK;KQ$h9~?cXD8yn)oy-Av7oxRL8BeXyl&77) zUfttMYCBZehEa$lv^kBtks+Y?3XK5tfFVV~Y;xR0=0`rC_8@9Z+4uTUU1&PL`~F(` zpHW5WO^Pu<9PiHiA6X(`Za9E0*A@sL4-}q+;LIMAzM@(lRbIpO#|5RIn=uM=Jv$Qp zvY~u~v)?k-meekGrW;*+P1qpD!0F+nK9;WsOO&nj5@XvW1>8e&7Q@4<13WV*o`u6- z4VFVD{!)Pbl}_{7?0g7CV!%{i5fcl8GWaxPCvC`*P>7^8h2RzVItbRV7DCfNF$`eo%8yv>z9l5{OL&__gcZUm z`rtYfC7)A|Hb-S1s!p^DS~?LVns`zn-LLdDA4TJFa?^c|V{O9?{jmJ;sMIrBQB_I~ zHNsY?_O}Kg{up&lxzWY4eWS9o$yo07&q(eUrn>rC|-#D6oLfLY6Y+gDvo?F)pVdeZ;%6Hx#)b*xwXcNydr3d$-CF zDr?aaB@>q2a_`)!z0QQe#3sD*tqn_yMgP9=*K117-_hEfqeFoa)lZ-O*+VnxnYF;#1~)3tpo9-eCpC3fC$=}gNZgGO%6^2V=dBy8PY?j4qGoc0-1dp zL6@iI_JrkE-bnr?02 zgVh%9#+NA1mq=8wc@~v66@O8cHEuQOcP5ozTelW@kDML9Bfd$(zpLqBhfG8~>~daO?}Nc5*0w>+4lJ+k;G z{)VTHi`}w0>-!6K8LyVTtkkoNevdlJ7`pZ;RDFMT+v(Gl_PcOrNy4_YS=rfFY`+$i zw7(?t-)-3)Z8#h)JpP&E+$xUN3yCJUzAh7bYcfi6Ilo~5EE>o0~f5Xv<%udyMom;2=L)HQT zk8jeBR$!}rMDwsoPTiU!N{~f9&SBLYbS!8l?EdGC{Jcd+TW5llj5*jd+i5&9Mr4Kqt|I(~W;Y z{Z`{ChvI3ZR*&>Bi}I3X?28Aqd7F;wCBAIE4)pg`uZV!og#W$3mV%H3iU26>`^)VM zYl-_w{*(q4tc;Ft{^@7c#o-o_Qw$+AG(r>cG}dsfXK5`1LfV8;0MHn+qUa2rPLWFf z6z~*W5`2092}hxRQ`rhC-5hC$K)=;#wR(QW|FGyd7=cnBuag}=f`hNeclAN98%rb| z{|5xty_Yp`;UVqt^Z`EOn-EUZr8${!A|h9@RUy0e`~Hxi?=&ep7k$P{zZt@!e*bvw z_)ww)A8L z_WB~Z%BUeIK1i_VAjA|qK$F5kX78nPf*He z&<<$Dq-W%`G38^Xj`2CKD-JS2v4t3R~!yZGJ{ zBj@!|13V?Pk3d0!*j@9Kzq$!^8>vR;no*i2TU=*ED6Zz(LTfI6u8hEN?=%vg%w#CD z%WE#AC^(lgq1+Y7p$~fV65-d#(ZSJa;c+j$v{=OZzOzczz$<8@Hw0d)LaS!q;Al;YI9Jd|9)R z>}*QIKj%XZH=HEuXxxKao$HK>(PQ>5{btuEFrwrSctP*^XHk;rJV06Bnoqu)7$X^h z2^wFDlh@XcBy|cDnC2w<- zGfd`x^9(2qsMC8sXTirHFAr~`h=6_1KXg+hP}9s{Q1zm5{=%a-{fmvh%GvL*vtntA z0&=MT*56zOueCM!#UD*5Jd!iyAs7kq(!IERlR`VPDTI_=B5LNbl5;7pfp4JpQsX6j z6yq&Oy*%vN295R}F3sN9bTzBJ;6!(3Lg$suw^FGQ#u~&M+xitIWwqz$1_3}%pbN9n0 z@7Fb;H<$!e;wB|QQ|SeXhEeOO;=aqSX1hSaQegR_L|Q8czg=W0>Uwqb?GV&L?+&45 z%V=l3w?6GA)qX66k4G?*mLJ-YJ3$bOKEB;=k&SIL3l}iYr`t1}G~|k34z(*nRQ~>? zQW0LbMFS()pd^TG2fL=Qs3UPp#AK+pox=EaPk=2F`luh_%mRyI= z*bay5a2}nP?~Tt_3M&W`nER?Zi)(cMdE~e-xDh~5L|^aM&b!Ig0~r;~Jym8YFJ@%5N;tr0AVJC&Iy^VM{&1Eao%W-o) z%XcQ5fAz@li_+~zAMZwEeuk`!R?Y&Io9jKG!Rf}LxMmv|ISUsC|1b=F6e2l+4-SoS z6>Sz&FPkE)sFk#*80E}USO|fO`QeltuI?_f#~@^3+|g}H3t?B$e$hzzVd{$XypjA; z)(-t!=@~aD_3yA!_4xBp`$yhmLkosaYdA&M)ee8j@RAyRCP5cOq9xkKI|WRg+PKTw z+VcH|w<#72<7q{SU7B#|9tIP_`hipf8_Sr2?VZnp{JCe25pZ-C{k%T-fVq8nAdD(q ziBMA{1$qHijP)MJ#%^b3t@EZ@}D24JZpMN>MoTzZX{ouhJ);=TYeNYpqkJiitMd~3mtf*8i@hPf2l zo3pA0J?*7sHlp3EfMqs9}&+W>}-6^`OgY8B| zF>I$Ba1;2ve(i$tO9U8IY;2MNLeQ?VwRHNim^}@+v|O<%WJSQa)RgWqMk%4SP}Z*+ zm+FOA%9z8ni%n#Dh!#2s~`>?+Azh9fjTZ(aAz2|%Z>uVzVB2y|v}zyRit-;ZCOwMfgt zUd!O^q<6b>vBS88kJNqBL?e%uv6n5N`zJule2o5-{Ju1EtaMc?`FQ*F({44H@Nf9c zj~?`hUlv#ae8#bus@OD6b#<8iXjuE~#Ll}r(!T%_#0=pUFQEuRGT&WLlriZ2x&pk* zE5gFTL8KA{04ocv&C9g`)(r460YU;A8-M`IS2Rt@5|vlzb*%i>>RiaOeKE%uC3nNq z%82{-#nDQ~b`bFhGx8YMHXhUh+LF)>MMa{KR>#up6%_&$gt71J_<=I{*2J{sQP-S>s|$iVqjNuNYJ_@Yiio8hb6WRVRPUnfI0v z5eeb<-|GUxc!imn3t*pC@^P2R#HM5q+Dw?ALTT|XECt-<$Nj33yiO10$Q4r}96_lB zf~W^gBI!y;BrA|WRsBe@^>QFks5~pO)94m}54(;6?opt`80kRx+_T7=G2 z&2U1?{j*%u)lKp%ROaB6%;BQiNx#g=N!@td)$)a>Mbjj)NZC-W=P#$_x0Ogzw#kLb zT3);>pQL!!69;gn&r-d0E@!L*8WYr|24rQk@tB(7-`x89+%Pr`I}Ll?j_ZIlpo@!> zw9FlTcTdskZj&Q|?M0Aed%xzp^10aY&(+-Gw-A2*rL#R6C|_f6?CI6GazbEB@06O>H_x%+sJ@wof~IFO@ws{V2eUbT@}b5 zd5r>Jkl1Kw13-MtAhE*7WQDZS5khbqjqU)A)q+LAV6T3K!+3=FgCAvVpbsa}4x~=! z-RV#dRA#_P5m^zpCO$}{jx?Wxv$qx{iIz-6HU7a$ zxBKvfK}tXWWtys#pWQ?)1DVW$fRKu97yZil@hv7h5m54O07)6387Vk8gh5{-0~9)X zg-8I~`&_rO4VTUk#0FFBb>LwuFX5R-Xv<*Lnl>qUp^Z5ZF6OI`{if%~nH1+B$ zvcaO*CZtKpp=34R?aT}vpq*KvXZQ6L1Ovj9Z=6kJnoswbXcC`K3qp_@S=Jy>sRsHK zD+4BpuX9!a#SOi#mdmiIhya7*1E>K+ay4qMfwmaYFn%9mOmz=%=Y3XXD=exTK5j<* zTwyuO5>Q2hoedQApcqqASh*$m5b8@CTm5RYPXs@M#6e|lwA8bS<+|bab2{Fa7;6-{iW>+9eMzNlP*iK^i zZI-ox)sO33@iq_H?_Gvn1<=Sq{+GX+y!F}5V7gE7TwAi`3;hkDPbgG3cm$O{ArJyR zBYs>V+DOmaplCAYTz?cC78Y0sT`yz&S?Js+=xS0)0la9q*i`=8n^oCMAYtF`y|%uw2)Z>F9k7WM&f;aw{b@{KSx%F)t)IDjPJu}YLf(De%~39nA& zGNBAt#FseEI23B7xoq*sY_vi*Q4Zb)T*k>#oLqP4HbuG#TCM3^K}_gXD8>u4qU#-< zpU`ce;aYheiAnaG5tpkm*+@E?H?p5M8Z$w3I1m~~FP3Qkb@L!)W#z|TQ#jbH?l>RI zEi+of{fohiy-pNLNl}gs>)Lh397FX*JUq;+4K~Ny{KK+!J;1_g#KO8ndB(BhCH*r! zE^c)k9?I~?L*IT+7SVsPqQ?VfdI>2&?RaRps~k^s3+nbyKY{?e;vO$P1uE5)9(D;p z|AW#}G`Gaz5mpcclxqs)SG1rL3D^DZQr-9yy8)tfXc)uf80R(T#Md{UsqRQW)WKnb zM)eLZY9JO4=`2cJIe=t%Zes5|`t%8(74z}hpM_V5TaoyU_i9K}&5Pvd^IHPq=d2Rr zW_KtI$yq%6vOrVD1CO&`2s=`0sm8eXY-Agth`$wJ``-v?=$84UQEN=G1!J76Mc8Tn zg~;wiL~#89Gd{dNpc!=O#bo%4=j9=6cX=MJGZHw1IdHP%HiZD#lN!KBT!1K!C`GHk zhc%&cF_3uXHJs#(iqt_=*3+mtP&FQ{01!dl&(r-vk%jf#4gt5WY(bR-Gh3BoJH`K< zM`9k>y&plyH&upXs%K=~k$%751aG1tCBa zP%u_O4d@`!x?$W*@oBEG^yXlofSlUFovcOUtZi(LAI5fKVMgORb$BW?leA*Nuq)kf zZ_WLS(A%YH@<4hA)+}mX(Sz$a@_~V*QM1o$?WNB8ZlMoWJsM`~AB9O!!i1n6B+dXK zk0nIp2Fmtk_(6HBmIBXbi&D=|-Z(I_p07*)K246f5qyv`mA?T54#-;}HQ#y-)&PAtc=i z1agL*bb_Y4;my72H%OA4uM@v39$OT+%+79;w~<+;>Cx3hw^8X&fAiIo=xFbjTu|)* z?9zbK%%vb?%VTtenf(xonaP8)+P^%)4x4cI2pjYu$tDWX+;x0qoszgcl|j>7QjlOY zmzd!p-1>wHlfwp|&}Lf*Ej18j#O>h*Gdk_TEaF$l7U&GvB^MJ}6iX!R**Va1Zg+;q zw2{`MqCktq4z9;W`=E>!q*Mp6V@l7o77)Jym3)m=WpaRyQij7fsPi|kG`537+$WT_ zc3tCj^+r%GT&oB8W|Mu(i_(SX`;tFci5e%_kJcJ|rR4CEJB+PR0@GI0?3=07^`R@WBMMo$=o4NV)p)*ePmge>N^N@omZNDuQg#L^ z@T7=%n$Koj#&ccG&=MZtMP7_o8XU`wFG`B!!S&Y(u%f$7503-D)R&Y;Q;&|1S|7xp z1D(Vv=^_y6Yfn#ZKK{ey`&q7qOiP=W@W$iu`mZ}>46P4nxA{_p)3zih_eiAye>d|o z^lbB)%W3C}*f*KohfQ|@d>sY_GIe8x0auu*@?xt#Ul$s$3jxg){ieAvDa38;`$|iN z?Q@VRLUwFy1R4&^&D+K*EO``g9yav&{Vsht`&4xOsiTj)v4R$eyV_hPeU%wP?#3v= z`6?hRhY`di84O^hf(m86(m9l2-IgP_iL?D{Xc4t=8GN#NJdog0)h^Z_(G0hz{_AKDt)a38IvSQ z4?ZNihR!+KkACScTm@aoH~U({-!(yO_HD`Y5mq4 zelT+x0Ug!6NHI=bYTRl=RpW6AAZ21~7Hd|u2l*=<;Z$YY;WycBrmLOtJOX7M1_k@M zuiXeX`f$jfR|dyT3^=T(_5+J2J*&3`Ej~3u1Zn!@$&(H;IjzUYs36%fATuiL6R-6- zh=S)vjsciKY3w7ulie0w7BcGo;1Crya{qIhcojrO8@Xz};IKB}TJTL)_(T%5Mc_#< zl`|{G;WbCZ5GS>;bzoA~q6))_SC4=5UF^ZubLLiG%RmWq<0 zwa3yI^HBpv2j57tti4vUtgq$=jC`)JhjuHX(fU&QEl`MDjTA3Y$^n2z&r&}@SGflL zm)2I-V~~^;egYZ6E(KxNKP;B77Bo}SMAtv~AEO-{@`htei@TnLjg7gcJq-+7y3nHx z=zd`RYHH@1JxzPs?99v&F>`2tBJ?J|`ra6zV}H*E`>0u=`yFV##-56m1m?1Y6`on5l8_7WZTNH z8l+@IOlP~>w4I#qY*ZF5(=#y<)%mjRsEogj!zM}3)Lx6G=bGusghULaH;F4jNXg^G z{wE%UlnDrh8tpyCuKWI#&H|2(g4^^&&-%Imz;GBx`NCv_vaF4}39fYifUX~c(irgc zSC<>}V2vfHZB`!VJsa(Pii6_Qvv4=F(epmQW5(33JRT->(&ttj(GVv+dI_`KSPu(6+dO{lo0tklx;S?F;E~HA( zx`rL^0ss>Pi20dhU5DFtuAYrftSwRsQ)cpSpFBrQ-lOkR>1Lz|#+^qrqmlEJuHID* z_#5xFw1w9XDPWIQxScy0e{rxoW^mIOu~YebW6*KaJzt^$lknRaKi2)mjPdhO7avEco=yn?-XG{=?GMMz`V0K@CS1iP5fY*7gk&D^S z(}EHd;paVUY%?YhBA!gRuMa&tgkwxS8s%uCwg-&S9Yyk@+%eW;i(rk0^5XVOtV_D^ z0?FBu{)K!4R(L!`sq*1dl0cnPbyYSJw*2aNNQ3)1FFD*AQ-irq^gO- z%;WK!h&7qu;owzj%uv|*+tqAVIe24}sq`Zk^4ad*L=Ax1Waeqms>I$6QF{ChTw1`t zgG1;%m}*IeitDXDPeQb~nsy(wv6%2laI#9u*Zkl~$a{_3c$QnX1WlsM+P1+l06keUzIggDT>`|;U5I|_Cj`HH zCxh#(U!Hc*xW^Aj4-s~MFqVL~yLf~46I z3KTSyeq)UKg?J3gI`ZrJB@~Ye5W2U`o*C1I2Z%V$d<}{@Cr_SCIgr|FnZN9X?roAm zvF5zuwUBilv*x`ROMvOKk~6(^<$Jro$&B;I`!c(ddq?eb$r4BV)W?aTZedud$G_!9 zKNcC0%T-y>DQk}@vu=osIjF~c1Np{aR?i<_x8A6HW_mc=aLgg<4K^=LU{Jv$y=KP0 zJ7tmxzLJAYA%Az@uXH*3@7;N!?1Lzaxo?A}XO0F2o0;}f9^{$`XSKtR0yU*XN*vbi zI(GdvWJgeNgPowgEmN7nnaTf7ud~)zk%@4R8a%g%_&T|}_wh1_lH(TTC%4Wvr5lCh zP_ejlT+oA0;j5ZD?MZ{~2YI9HkGM>TtrNKaAu5E3@?k(Qs3Ib6Q6N63Uw=y<*(=d9 zHs)}Y<}ijbyAVJLBUH_AGY+ZL%nA?^Ua}EI^cvyKZ;gl|zGjwhPmEK*5)P&_AzO$V z;)RJg*#v~{uWe-2o%nBIKDS+RNA^38=V{}`e&b7+KmXaYU4@hVlE1rrZb1kM7*rO+ z)jy{GwTVM#fQpIIrlshm!YVmq5-j9CL3CEnXQJ)op9Tlv3c-JNp#zmL*-OAr1W_|A zIz3hp0?W9zL3{ua{R_(xM6tfjAc-P?>C1+w=z%pB@$P9tzi==GEfJ(cT=|g4I_h9` z%u`16ui@a3fnQco!X-njx9Y8=5Su(O^*W6>AlV{DcI9AL>%^?+uu2H^-HMe!k5vJ3CIe1wp-pU?^rmN)ECjG#7}McK`fEmIAmA+{iIe#S1qtQ=Wl8fyZ=73s1=T z4mm89eqhEzyXr+8fwbN#50zOGDBIECG6x1;JR+^DQfp5)x^A~jb_w_c>(>r@ZmpB$ zBNF?UC?|^lMyVioSi587w>oo$DjWy0S#V@m!m|rA2>v%md{q!#iq42gv2ZFP=d0lN zFhDq?R4;avlea`8MqJ~v26~eqp?9kCD8P`u4@7&aUc>e&P&amn$AW+FonF(VOP>-t z39UU?5`g!Z*)gv3H61;9yR&n{i6ydkKmas+fIEMhNUq&?V`T0)(hk@80;pxe0lnrt zMf7#oB>$l1?+~I^Z=K@Gxt^|?4ayp7;2D|o2Etxyg^iHwARSY#-7M4EXySx6$AKev z%@1tS7wm>^k0mw)@dr|FPc$5F<4dX5p=4hU?DJK`o~!)Q6B3rAi)=z*1KehPE@k}7 z{X~ZJJ83txwb!_A=W>fD+;`r^fmIcZ;h3QDc?VS!j52O=2PwHQ{=wqm%(2`~YS)|> zrlo*GsauW&n!ug0WMY68I^vonbyd#u$ssYSpo^^ zctCl&RHx}`Hkk;anfvuj&Kp+m)EeqsqChu*YjTG+;Q2~;_>seLK*?-4;O4iTUv%=s ztd{%`57RhPH-%)*8mp%k0Fycb1LcjzyHEV}2(d{AI>olx%89_WQrh$z$%$Bv1xN&= z?}Xipfx-|^m*kCb26CmoaZc6?gSP8KV&sD21o5tEZ7YYPc)0Rh)kQ9h&gbL~cDN&~ z{tSih=y+sfZ#S%&Eui$|x-FzEAzuN@yjV3o_n-z!Fiu71bFL0N!Tbi?yTT<(u6m{;|(Jk|lvih@e6;E_g;&&RUQ24^^B{To9^P?*Lyr}4Cp|S!Nh8g33>liE4nrB;r+`M6l?F2248*I zd6vUKDmH~R<~g5^z*J!GQ~`(5?OK~460BOAaZHj@E;>O|=eo{J>h%?`{QIxLYz%clOIrcRSgDa=81>EPx`@&~YNd742MMak z2XC1RbtIFX{3YrTL~1b-;V)^;Qut5KBIv&TwOBi;42nV=fD5`+19^(!1E`61Tp6S?@r-M=(Io2M&Qi&U#=AKonaMo zw(@Uj2C9hv?7PsKz;OVahtiu&>y+9QAK2_p`c%m|ui&%-S~Td!5LSFYUo+ zaU{Zk6yEDBJRS-e6D;VxSib?*Q`xpj6;T0)oKs4UvK<=lKfWN;n-Uw82>yuD^5!CRyxF7|f#X zV}B5Y5V&@Y7z8)_5hLTAr8>AnskrRXFqVV0(um zt(rGo;L`%ecdtak%e({{+&*zYXpW7n?TskWaUmA1Nb>^lhhYj7Ni zMa_19Hyd?VS$oxC14FjSrPjN=VY|JrwN`%@G>hpO7@A-rJprSLaa)=HtMK{$1TqSE z1~P<_iPVlxN*Q3?Qdp9bbp}8b%P^3j$eqFtFAcKhK1;gjcGDCKw0C}()&byB;oDAR znHQ?+-VMs$|FN=G$kwCN@&S9+xEUbEnEQ$Li$19hMK~}17a@E z3{@7a0CHq583_eOklw4&;!K@aoe-RYp(u zj5Cus)rDkby!sFq)Q4~%f78ITP!*vWuOHaI+Hc(7w+AF#6bE4MhA{89J2{KKcuepk*??h@oEDNr zX*-V(k2CER6sr?Vv1X*=9P>GprffMBiM=_p_183W?YTWOl-@48^YkeG+dsC2+PWRq z@maXgT$$>&*-vcljk=>PMF5)WralU5^0bERr3EF4G1er41)lw7Ap#8V&-XgZ+MXNq ztVz?3iIek6dI{QnD}fW)#VLe?wQARS;d9)ef$;ZP=5zB5p1ZTOLf2s7v2p>?9kuSc zHN&Nvi5O1W!wl=(H_Q+THHJc6EGGTDSp|MyXSuy(S(Tak{FvI9B3nML*b+3}@}S@RQ437*Z7rMUL|47bRMtE(_*AE|yYfSWC{dR%|@HCL}NK zTGcY}%VHs1Y9+ww-ks?9{U78lK~ zvDPo?@{}Y$O&m7^IUK9U0MDnx5|EN2X%VPe!zBDYCT~)@murLIva7PRaKgfCs3Nrp zbw5l0=uXagaI~hZ7s_Wfn%QvJ{55>dq?g9`N}N(PX(7F80HrWH7mWz>bppmE_3nKw zzEn0lquiY{_4Asmb$K6#>VKCf;4@6Y#Mb9p9x5u@FJ;d|4tLm;cYs)LHLqN@7m9`| zwOd$Db?g4Qfm3fmI=*Nfqd=wXf0s>nFZ7h5V?u(KUh5`JhR78 z1JtqC=H*}hccObOdch6(x!`yC7Z4E1W=5zMfdsi53v54S4+d|e1hQx$#v*;;wZsV^ zl|d%F^Aj?Vg?5^MPaSA0DToNX@g^ej25W3@kB-mkbqHCrtl8GoILhDNj29(-Wp~+t zUFP!c#H=Bou3Pq~GxHjkU+%odl5a%&EN+H0HBDpf!P>k$PhhpDWCcsJKk9Yc z6RB65d!Z}+y+gaCN&?`cdF_T+RHMy_ zMi;N%fa3*jY%2c11>&x$s>d3zaZDXK&6)GInb##H6kY=Ai1fqEZAPG&l+wE(B3xOzX zsnV2syLJM~(vUCJ)o)pZSl zl+~!{7WMEai&EouC>uNK5b8r$Lppe?31Dl^tR-d8 z-*EobEqV`-Klmr|zztx%bRt1=chta-Y|v%w0`J^zWM&bbiE3~lYYz_fDd z8|7dnpsj24S9ci~ORvm|0uUhJvzh~rtS40_ufYD+8Sg`OkO~tZ&dFnx@VQRTfII&C zC~|)E$@xbeTOD&QD%HmVDOuAd$J-{Sus}*OgE(`7E-?r2`mz98ud6c1ELUg zi`~*ERB0mHxEP~YEev8T-50^+wdSX7ejO0M_HlT8C*02^t;UC8Je64H8{ zg9z;+L82}+`N4|;-Ze0vGM&y3ekkE7x>(3T=xUS*0cXe^SE|v9+oWWR7DwmNxy|vY z?td|{8!yx{F2vq^=r6hYlSB(738$(ACG`^KXddw>K)znLh8J~*xpZt9AU*%xGcR-} z(jl4Gr>=ph!mB?dFAr3zwb6urP_r&vHllFqb*sY+p-|;WpANoi5|HVfJ=c)yH6sNo zbmZG-B&3Gn*auf$h*3!Z5pxos3qNkqkC&jja=g-WkN<4cGmyD{BkTfVGVhpzDr8sX%xL?L8Q_EbGms7n(iG z=}AeCFdS8oDM2pFvWTc1In?i2{H)(!H6GQsOb?EYjSXBGVWSLRXHGhCXZ;_O9@eFy za~q9|pfNa8A%-+{r1;8ut_MUia*F28=!^XKq999J}9McZ`J(pw^z6X z(67GTJf;Idjqr%*KH8up4%8zni3mcstTbU1*YhX6uO66}Rx z6aG8fbI8I|hzU9zAn1iU7XdDsEXIEl^%mOB$gIs!dZQsL>y8b*3@>c|_LG7D{QY)r z@lU)KfuupWj@#elGABAbGY}~|8s5LyF=cH`k#mparr)wQ)Jz&7y@^WU8L zF;wSL2B$8E7%tRA>G6D!g=fY6S9iej;QFKlRmYRJ#$Ci^hS+K3R04($aDctwjWd@D zr|X4eFhdOQNs| z?)}>j?uP@4#{wIHdGTX4J%T&=?FLRGz=1m za^1oq0BlGSqh?w78;zr%1TYRN|2Cb$1)wr5jRuG10F(~d*TQVJmhy`b{g%b|z)v!c zog`$%ocfgK00an114!{fTgukl9OrM^8xQ7NfN8{n zyy(Ab*1r||=k2d?$0A)Eray0xf2kaW2VQA4MU0e?um890P8f~_m>;)AQgN{m$rkr^ z=kKW+bfYEriQsdg#K9`__epa2PLel3ws)0anXKTwgC4`=lbR&3n3>qApXUWIYDMiC zifwSHMxHyzA_DWA7sr}D1Szgb{R`sr0LjRaTn*jyId&K}=fA^-Hu`6xoFnt^4Y?6@ zIRN6E0`32AiY(DSA8#xG)acQ_ja>L~+W-6u8N{9(g*trt_dVIr$Nu|`e+R|jCuimp zJVX@a7|=+8#=>#?tTU8{M=R!1DJjro|7KX8g?<8^I$#~PY|25;u`&{RKJ1kkU* zy$D5`H86+rBL6$NwXVdLpWxC#GFe_>MKyijE{yKY{9H?#}o&1(AKM0T} z{yl>>K-{2C1r%H#X%lTDJaoGp4omAgvatS5D6~F63BlD12m>0No^EZZUO98|GS;a^ z9+*a7|J|+`Wwog6_6dyE9>e%L?E_|5X!iP$j3gP4N=*&DWW8ZD<0+WWbU>As~t+|=YNS2pJ%VTPv z4wg4qvQU_?0{-o4rZV^QKTmRmK>%3=aeRWa6Tt@~%S!>u)VkRhltHCJVm0ADE^Pqy z$8l(Bi_JaPgii&yB1h40n1C?*>_oVWmc(a-BX)*`z~xDKKzQ4V8-8^2(igf6bT(y1 z6A&LlAtEFLH6gvu$;)WL|2_l;w+eF>MeI*5htncZo4-*mcDU8QG4`Nm8e|Dcz~*l; z(2cd@F8z%KrzghuMl0xEgFij=o>-*puDdz~g-oa~BSwz92lFx}z3@In#z0mT5MK^* zl4!wy>rYtoH3OG&e}TBKG;O5mhbQUT**({O5l?Z4=&8kkuWqM`k=!Q{#I#G9UDACw zH0f--T*&Xp_SuJ!cddh>Txb6JXUVlG$Q*bjzlZ0b|Cyj|KP$H#0;)Dr)w@Ub1-7L3 z;8#?*4Eda0mh5G?QEe8rPS4pB7u|Qzd+<5ecd6hYc6l;@glO#!5FG69-UbBnSu0i%W_czU(u@7 zQ&N$zovX0iKuNN{EGlR4d>~033jO?rm?C5z{^9T5U8xo)?}^N%fiComlY^){ZdQz~ z=*S%dy2msiEN3*>u6hO6;;kTC*ryOgWf1`CV6eK${=%Y%>&mrjJfgUNRvi5!8;vP{Y%p|K;-qv8TU5OJEuk$@eT3TRXKJrA2gN@XCr| z&CAod=4N_=eIr>UIOp0cH*pMvDn~fj5zzA^`65a|uQ0}QP|5~sDm2YheI*go6c7}K z3A6IC=jMRzLhu|ySe}~7JfMt3bPUw#tukq$@sAuS5s*!y+jTxkqi;zCxah_LTo{x9 zY8ycP&$n^t+iGHF5W(hdsZwiktcLY5t6(bvIa|D+TGrgyC~M+-BmE{wAOpdg_M2q;LA z-b3+#7Jm2r$b^o~Qes`_A*Ewh9kYw+jch5WX%*-=0c~_nJ9AHrT zwynAXR*CBbbUaNp*_kB&(^Y{;v#bVmY?#%{>prqzf>=DH%CpsBdw;R}esuG)z7y2S zQKSv;5gu#CTs1TVXhZK42}ZsWip*d1cv7G#Ri~A|;fASz!|s2t@p7ay;2kznc=)zG z^`HF7&`QQxSt!|_>4ug*la$|j>M6zKm)!?lQHeh|pRUuoU`BT|w0%mbKjPhaO<`PM zVa&Fj3Hm~Z=_$XWm&HQ4(O!tRf};2P^5shaHckD_UjWHUZxjxWI=x7uRQH<|$59Up z^{coq=&b=>TD;vREJkHtneKhYi^C7(;K@9jfX;7aOTbU!_ygUW0PBX^9NokBnIIgn zQ>V=p`tw;qaZt4PGdio24uB5u%jZ0(Q2jJ(O6MbR7Q>E$e_RgYI^R3DA!niP{-{SR zRE*3#BbRbIsSA*FUj{`!N!+oL$w8G>1CBG1-Nw?sAV_gb_f2xnwlD&&jTHKN)MOctvg&F zqhREmd|17?hOoP=KwUlkU+3LWyA8PkwQBU9{RaoP`g(h6uex>Yjop$Gxgsv;b{hRK zXYdi{;p^%Nidn;_6dcWiLP9jsTb>6`9J+ z*BbQj)oJYnyLjE{j}0Qh7X+PhwRgi1tbIqaaf-JNvi^M0tMVCjf6d*H@QbqS+4qRR z`(;tFV9dSai@`KOSuB!X+XoG?C4VJ&9{h6p^dE0g?^v0{8;Vr9rV!9gB_4MJBFZLg z9zbQ-i?AIC>Y+*o1Onmoo#QmURq)U0m*`-}Qm|KJmvU3;-pKR{(<5s`#!Wuq8}&^a zfu+BXX1)ar-!SJ$SaN&-)3aUm5lg%^}+^v(u( z{h;*GJiXhZdigUkla52shK4#cs)zGgkwR4LQOnu^Uxj5o%B}vpD$(}}BC6PP5>Ibf zGjLZ#-|)TS&>6k`Lc7>`xYcP0ZuKc4`Ip~CS)D5UE@$nK#YMZ_Nq-Ohm4;oc15fbq z1c(o}Hw{6PoVt!Qqan=nBb2P8(9-6u%|@?s1A1)i8B^iAebO6<$`_mJSvyMLUv{2K zhW-uL#yY!qL|&TJ89u(I;?-l34eAyvzaxHBqgI;c(g~QG>gvvN$hkCozQh*c+f|-+ zSx?zVI)u-edG6-IKC!IE@^wBjPiSI5R5Z#_`2+5ZF0aqwaEp!)`W`8ZY*BF-ZDm?l z<_{s|<)cN@O1$BE?OVA1M3y?@R>QaeQ9=%>`EzdX)BUgz6D%;e`AN)!MlC{L1cA2J z_hT0?<_gcsh}wf%N35JSxGYVjqsP1lIw*fK9`y>ut^6|jL%^-+{UGq#gi^gfu)#?b z;D;fyBsuLK(mw`WO#(`1#?V%wy7aFWO_%Y87?rLuh1%>s;W^1m8{qu|k7()FUIfC- zf-W2i1JPlsY;}y8y@D){7vf!)Ok%y;&Oab4O^2@c&ReZ+ZY`=<6bc?~lsiMbF{kLF zT>Pf7sA#l=IUkY_Tjemyn%7>rvFHN)o zCSFZHrYf-0b|J?3prNy@+u-bNnm;??ZlW?L&)T2bKk<54Z5pKyb2?fe0V7*Vt)5G| zTrPrJOTlqADlg}1`nKnR6ZHxh$`K3|Q;JZ35e^BW*(zM1_(1Zmi)p=1(oud#p1%I0 z#=G_wkgyVz;Sp#%9{rwa{t@G{ijuz^vrlQmuZf7;=vn?zP`>)PbPkrmt7ksr0>nig zk4FXnbVRM?gz~?lOaHdo(IPzn{ZQnym`Wc9G50P-Q6Rd^qp>L2_<4IC$IV!3s7)k< zs*~eh47e7`kR{KjL2ao~EYv~Dw1F?|3EvzuGwJ7(mzQ^#JhTWcauO6_Qpdm_#=#d{ z&To0wzwtWKcn&f=ax<@DSzO1kvch5wi1fGryTLleNG71prXDrD_$K$F|6gwLZmEBTnl zJzNW_--Hj;%|D)YkUZps=r(q=0TOJOn8-7x^;ed=JYU^oLp+g8BJ>hw{t?@6n^#P2 zRQR)mTYusDAcftcOov`5wF*wY5acOvtc|_r?D0SCkPga++Z|jRrSf=w%Z3lQj$lAi zEOD$(u#djWbE_#Y)CJdN(n`Es+6$@~EjA91Dh<5i9ZSX8k$DT|2tFAA)5GS1P&G+zkewz&ONE;C|(~i z`8972N`_^kA6xafRx>HuPM*T+TSn}Y?yoxCiI`x=*wemK%d>MmtChTZPc@aUx;l5O z#ch0`Z{fN6D3+n}<;M@5V~e+XT;98IZyK%2Y|W>E>7e&RGAnu@pq<}tyoW2=Ej0Pd zZO0_xfFf3*Xt&nsWr2W%!d~%feU~dZ6n$MjRRdICwn&o5wYqt;RwPr{PIw&nn7v%NDQL{dPK23L z4F?!jp4BQGvu#W#eQmB>Nbp7F;ZV07qb#5H?%E6y{F(=qyOMn=I9p>FTE|ss~ zFc!4_r^_#L!q+-L<7Q59FwC3wJ-#`GTQ+51y9c)V!}=p$IJIx(#h=;lAlpn%3Cz0{ z?RJCBk4@tgsax>mbaYsyWpllPqM{QM2@*+yh1bKAt~8$jyV;*R`PKKz>1Z#QSkJC` zjHp-(6qopRi-%T49g0<$zAY!O5ux(1zf5QTRhH&2MvIweMFcV5y1Vqn`y@J;r48n* zu}b$D6U{fc&pJOK;q$6MvIw#0q8dKgZPZ*zvuJ6#d3u?;lLJi7Mc<4jRC3-epE8R! z7xh?bgW{K6N6RcSNA;easzq2Tz@=k!hqgA?U4HbqJRZ`-$Zo4y`o6-=luU=9p!I;+ z%AT!ZKw0l9Uu280{{&g<|2oOvdD24i!0b+jO)Tnbo?GRl7OHX3U=g|SvQj4lQ+^ZHlbwJnoR~%(a|KQem zhgh)7Fh;ecuN(P#ABL{*;rXf7hGt?$n*!F?j0Choi{eE)#5StLo`;E+!ZJTP1(GT{ zDbo2-Y*9Jgub(WkFYU=1&Jn$5vZw?di>;-mw0W!=>i#ocX6i3U%>KIn`E&HM9S9-o zz>uwlIGW5V^c|I4#qq3)(XV@U`IonUsUWue7lz9HU>O96fbAuOlFJut=?~|v9%e%2 z3Y*3p*Va_5pJVCEAC9PXlcLe&sH!^l*rkq(4S;~N92L0UTrvjjKJQ0k+X&g^{wKh- zIPG%gTlt6al?RvS{HRMZ9T5dxP*zu<|FEO?jQa z#c{mZDT|`_F}WI$@RmM^h|dBiC;xm1I^@aw=RB@D5_>Y!yTS3#Y7-7C>v(!RNOKmt zAnGbF#IyK?i2>EvM1f*vi#%Wy%pk>aa9`)HisN^X2PRrUeAfC$2w4#}X8_aBCkQ)O zER{TCy2^Ve+O5Gh#?rU%*QHv;Wj`G$wjKWSewMpjX^+<_ec|-MH!eGw#~(Jj#h#`& ztJnsrunnSibJ2_IUFTp!x71X=NKu-!2e4&enk9Sy@S5w3ebSZ@o`;>XU>t@-NQ)g( zrl$iVW9PO<=}h<@F85dlR^){&-K#Z-t3BwpG$yW71>l$cUY& z62fvl7CQg`9Jg4OFM=@15Hyg@DTR>*5va8!R45sy7t9UI5R9jw*{pZZxJ5z7#a%(u z;G)dc;x=mMVB)>ysb<)(q7gBBxN?uFU5PrkQt1oNTI?Sdh35+jp(jAcH(dB`8vgL< zo|tfI5r_5g8tzJ4?tTfFk~6G;B^=t2z|s_FU+J;UzPK<{?ccPzK{udpsY=t{bM)X~ z>~bl!PBL~j&<{yjvqSB=-3N0Y?N{m3f?#g(71^|qoHSQ$4hM%*$T^XH5s zz1&B6J?4e{Fo4q|k4A@q`t$eqXYAXHFcVqm%2FmYy_m)>mo#9?&W~J-V**^x19}$r zUXw0e#zc3Qh-GIm^Tn;Tr~JGOBSVVtzI2RSsOU+u1(U+AWgHF$KMITWR9#meSzd2Q zA&+aFsBPG4N?@9h{00SR_^S3BD}-^!0#M>JY~ z_1S9pGDm#z1^^$yP4}EVKlM%fhXtLV=zk~?`EqxRz&iI?xYuYw*guHq)AIZ|dAGoq zNcLCV?KpWU?3Js!5a5coOMok$xFz&WzZTa7B>^i%>b{_5P2Ri=?-Jk8Jx}Try0O)1 zcwT9lbv;Ve9>-KYa3LD3ldGT)v^_hfyBpCjQ3772S(%R^-3#>oX%Kt3ajJcM5qu3& zBVT%+j>>^lk!xC^cb&VlPd;0kg}h1Y@S*!BUgzCv{CzCb#sj+Z{S1Q!Blp%O3;INj zGA&}RF@oT;IKWeJ)nsO2p(>->!uPO-^0QKZnW31_xdiylPt3QM={Sg0D1YSo74C~O zD#T|@i7gNn2m|$=d^kS;lo)uxjhyr)T8c)6Ss^)JiE9Q01r4`;hW9Wgl9zTJvoUTw zuX19DZkuH3{#4eU6uC7RbHGTPPA&Q<#Z~NZg~R$2SV{b`pK?NtMM;Tfur0rWazP-3 zKIQEQ6u*-`>Gb>1*nxd3?z;mNn*f??wcXTAY`N6nmn|X?WG!BU94&dqK>K{GvFKj= z0qgpz2XI?KzAeHC#J&?=$bk8eBG--w*Gh1{!Q(bR`n%}@SV~PDE+_lZKXNgu>~q(l zTTHO%Db^B^fGeDK9vpAafV=2yZ+ynI?M8r=mL2~NIw(A8VI4EN=tBpkL4!qlgg%#= zb9)i}*}}x^aJ%JpYbRC4(cnZGVUTg&J7D2DO`^{5K3cs{B!IAXKjyIN6p} zVf9Tw?>6&fd3Ei7FP=Plh!Mo%qW?b=2bR`P>$@kxiEh`zesmxcKCcj_;PWjh zyZ7Sq?w!or0=d_kM0RFf^N#fogWHAwKW0{J>3fUS-0L7nIviP%*I1oEdv3!+SUCr& z?LWH)H0oqoNk`>jNg@1NZMudkB@@x$IlyCoJ@IQd7i~0-rWHH$T&NxsyYNZKuZ~o8 zh3=e#w8y#~aXQ564KW{?W$)=neO!XZQ6wizJF)Gz^T*-a&V>UEtLc4f>3&UGN%q|v zeI_L{EzAO&EZZL$Aw{KAW-HJ9m^TAS8g=(MJLv6L`_?$T4swdN*Nnqq_ygd-5c$rx z$zN_73h?oPy9-h(5Woka2^|y*$)s*R`u>j{{mB+=Tok^%G0tC8-rpI|On`#lvM23ptT937k6 zIj`Cn2e>(lHphG8byIpnu^;=UoA$6ryB)3-kI~s#vBw?!TF8kH++BrPIA`f8w5qy8 zA<%yLJzZw)Cipx1qfJN7RdVFyS*lUyN>-jeyU+=k8BliXVOTpr>tNN>rFI|+A1cpB zHvNE{PF0iUiTif#1JoG_OQ6_fgXC)rAD`BMH^+H6PpznCa(yvYFq8XS#m^MD0!>t7H0IW5G?Viniw+k_iYHxY4Pp zaEZE>mXPoFOJ{DBPv06Mvv7hZP3rc$=o_cIOKCUf_PtTEb1Z!sjwiQ$0avMQ^k#za zVQTa{o8y8iA+kXN#Oz@|6I3qDP7RxqN)P!K)MCe9^C3w`)ZNGQ0;of{k-=``EY)Wb zW?Ss5b$5R!$sgYwDsOO=Jt$<(0e)(vhS70w?zBSrqYu(QUh!%|q_cS`)G;t}-HXfJ z*v9y3ER}J-*x`Vk515^XC!m$P*-70+SfNX_2i*=%cBsQP#-d*@+v77l+*>Sv0-4A5 zRl9e|j~(pW=W~YhWoa9MVh$9CC6lG{6D73T+^3S15Bfj|M8_aRp4T2P-FLqc3^2Il z+JZ;w=1cdpeZQ~Bey}-LpRLaF!dt~Ac8j$^+3?wiBGm&~qiTEyuh`)$Eg|Lqsg5EK zg1%P0Z`bj)W$9fC^ot^XC5c~T;a#W zv6NZp1&ZHb(`dU{IsN3un7%Z?-+a5t-J2z7?1@@FKCti_>iCoXq6UNlCuG2N<8Zh` zek}>34DVRtz`($zEKumM#$nxQmTI7-EpqjI%bV8;W9v`%V8mGzg6gxi>EL?pQCV*l z@ezm9mA0GU0VfFK+3Bn4TE)U)XK;Xtf5E4$Xxb5v)PX3{l`>Ikw04StxJOFJ{S$Ex z(ESh&(abV?P1f4kZYj>*)kZJT$24$h1-0`p`tIP=o^d?#^5qF`ZY?$q&^1(L`8zcd zYlk!4kuQ-tQG&;?t`%|xS+Vf5_2qKq{*lJZ*=1ddy5B6my$k!1lNxkK|2VV1 z%eaNh(2RhpOkBa7`QX#exMhUJmDBn-a%Q@oGf0Q>+vD%p?R}5E&S&PI0WkTqmfH=^ zmFDtuKeL_Cm$pT&LV7|G9J1x8jaC{YL<;7Xaa+r7aWKEV%SwVZlknm#12^)`%O3*+ z@mFfEhX5`B9QbA3qdA##EjzP&?x^Q*VY<~tskf3m!DQw!`p<_Trvej%%xqwi;l@vA zg%__BGYp^Z#H^B9xN_$;hoUJfNqs0yHVwC@Oa{vxmg82En6NtaK=FpR?2!m4=!(YC&IURQ`Q*0<$Y3T#(n*hEu0hgGWnemSuU4PT_AVWWoLI8;V7e!Rj z1C0z6L6k*#^#Fzi)O8G>yZt<}zQjl+ZofjX zANgODZ}Z65#$Ow+bI9q*m7e2Cd7g=^qT8ReQ`;M|{lYSBJ1+yT^8@JS435sTbejW* zW!ir{%(Kt-!Rr+G?mmheJz%16$&cGY;Cg~9zXp3aP7C- zj@5S|)fFZRxzxXQr}LtVbW}stD_hlFj1YN&MFQQ-$JaKByy7jfjn+O{4Zy_4*puPZ_L(d00g@ZJ!S>IW-CSS>tOIX5 znkpnlYYXaJ+}?MH=z2&`Kwyu$`h_IgY~!B?QO-N*O7BlX1QL#j7Tqy2`tUq?^0 za-23ajOl)jBlNKYBEKBctuMf>wtb!~ZiVzN%fOHX;k)uNKzREdt+0Ox^SWAYL$6NI;6-#wfyYtsBJL}a$?<*mTdyut(FyBsARtS_HvJQk0H58q` zx}t;>-2VQE=6`mLo#tP-^WWSM|5sZp^)CeS0)*-vd(F^r;H&lh#=3`!u~D869|~)Q z*6T2~F){1{>?-{qU|&G=rBM!?rmuO2HS}#ys*&StA_FxO|J_fw z+j^0SVErH&^yi3V4YG0<)l9Y|zH;Brzh~PPUDzxwcJ5~@8R!TM^c_7sART?44{B7^ zA?zS~e0xEl%!cmw7^L+h?N6SXs8%(0WM6x*{T3pM;++Nb7O$nmEFI2zyR@$B;X_~c z{sjevBA%RJfEViVq8!%zT89B~z~h*&85GC-by%_l6d1wA{KMgvEOh|+@p`z8q|W?( zlTD57bGvp;IE@f3G7Bzzx4|d>jt)yr4KLq-`hML|)arXw2;IC5{MobZ6Q`zz3J2i3 zkJ=tz*#2%mA4MxG9vtX4WVE7?*5NY&e9+S#UIdBj$Y_<{9kSFWSRMl z$yQ3|3$RXG(Eu<1V*+E(#Dncqr`10k;0E+)r)>+AAsa;um;J;|J zXy)<$eHfK5bmci~jHM?3u0g`d|EbT_UPGXQ4?-2L$F|?R{S?CQ^j8(P^y$7W9RGbx zzx}oUvv*Ji6Jr^LVg*3>v_KR;S=*cDJKbN6g8cSFMM#qZJ*3fj!XS407G#5B$WGuv zf`jubW!&$P!P>J-qE1_`_khYZx)uf3($tXx@?QH5{Ppcb$*f7=wmq}tohZErh%f@g zfUwv1DPYEk;ZLAHUcS^$X^ZJNCf&Qn`MH9BANTq_(3IdyL*7l9 z!i@Y8i&T|$<^kul`mafw$ctg%`kBo}_nfMYQ+^kUZEXK?OALtbCvb(iYQKR~A+ve- zolVz5BqRfpM2LQ?u&eRO^oT;!e9`I+nW*MfdwhO!mtaFvwj%7HR`~F4%=U>oUqHB0 zGf}_Iy^5pFJ)m$rqF1qLXk_FHi1n!ZS*_43_+*_5RaB>cwN==b2=SQ3R2%X=!11caO&C8bcSUuF zn&XR~+0l6P+Ujl z$lz?@mp21#7lZn8`MzTQio5iLVV0`nE|Ew2TeG}d_g9ZUm5)qjv6_n>Kx*Z0;4WVY*wr}zppN;{j^7zQ-4&&AF4SFy|PRyw~nG)Y1 z0DPSK;MT{6h5|Jl;1hrg6v6P<3c`Na!!@Uf`rPW0uGA62)$PCqu(WADIp?aiw!y9n z6ADtVnG;d-p(vS0{;HqKc_1njT_u;D-4~BECtqQOj5LhycGNk*bZMo{Y@8WO!Lrw_NN&Flg@u{gNhw2K8ch!mz)k|49++*9t2Ehb*WN|Q8js{=X2hE)A z+a4q7(61gbm`R?#ePc^E^i$62Gi-s#HM%o#`R$COwjDy%QP8e!pN;X z6g6LQ)!zM)wz;I5j^L4Rq@V8c=;ZBZ=z19TUH?R3v7i_MgRBC; zd8u=ku=GSQ9QjZ2Oc77qGDc6{0H0@cfxP|KPU&Y%5iipeA`~GfS!wC&fU(*ByQs{6 zfbq!xqh^edhF}q+Q2IY6T)I{JzZ%3g7Fkbc5{DJvb(B+C)@@=Ozk6p92WdzWhLiZ6 z@8{~DrCH8a&rIemh9|E-%0NZ@!1p%D=D3MV)%<|?S(u1yqgOcSKTtxfAyB*0IQ6Md zvIhF-4>&nh$O%Cx4T5gM0rJG0?U*P~kWy3W2(0eadmE@tzw3+_r-y|9;b<$?rAGhp z?p3T%*XIWZ2r|uV8sFaDDyXERXYl+FS!>U?AjcK$o)e9*k@t)}f6w0qjY}7K5n@?A z?~)$4ywHAi2i{bUUyYlVuKy>Fobo>YB<{Jr{SD{^x8&qwFoR5>v-jeNdwwlUus~^ zunKu)y@M{OI4ZS6(A{(?$!W*Lm#*X2LhlOlNm|!`O-*%OicS!|$0uG?TS{;eK$Dd9sF@AP1L$!vPJZPe`)#YNQ_qzg1GA90D=dOM}^ z-E+v8*Ws~`B(u(VM9|C!(q3{#MR#g8S@Fbje90SS1iIS(4+}fxerCryU#X`bb%Q>s zf!m1e9c|O~#UkZ>hKdN;wOf$nZGScNu2s|`t=;2|NI-&ATpB5_37&svZ@&Hdr6rNh z3-eQpC51drVRF<;ON6bp5JuW=@qW5VuV&ikY158UV!a9F=AQDA4M8@Vv=?_KS_j&i zq}OJOYliag7s||io%$YT28SWO#!s4bc`6m|1I~zzvFz$iyBxMHpR^`6`Dx^bV(%V7 zxM-o+4>@IK*bt8C2N`^eM)lcV<#~#&YszNYV0ft%->i` z4aYO(eDq0oNSpjv<#Bg0lij@5h;q=TTdIdA-JOWqP~@(k+Pp=uL}>nZ0n~jIM3_cb9ojG!JX9!aZ)X)AzP zGmAgfi}#GBUypf4X!J0MKYntXnw%sfQ+8ZHOQJy3oY1Ry&x2=k@1cgA-uAqM3yCrd zHlu;!j50KM*64PLC8BnL>3qt3-o}PEmbN+M5bMr+u1BSqu(0Eyq4%-H>gWvA{6}xR zx6M?lxEo(ULeXQkMc>bj%fT#!gBy?SdXX`)1~=li2C(D@Esobe)Ebv@ z8*jel6h)C)r!w;g4e}-plBikEVa{$8nnAJBRbsaa^M=P-?yofYx7K6h1sMc1uV1ZM z6klP`nD?NV(Q@W+-pIpK`sTI$?FvC*^8JaFv*(w^l~}4DbY!=np?3kv+JGEC-oA+c zh)c33dVjUgIa>MR)6}wMhmfK!2EY8)**A{tk@w!_N*H<^^*u94d3o{ru%+H`zXwDd{9Ko07%iv~jIuUG`+hv=NuGt)0nM>FoRf%6bFs zL|HCI1+2v$KVKAHFlmH~p%rEFs_+hoTbjf^#SCDl*1neTdX#}01W*e#=|?R(4^Ok_ z=QHN@c5c`kOvg9H338(qHr;(r&)WBH_%r5cW}Qmoox;$(kkbQd*9x65FI(VqNi~G{ zuTi?>@$s{aIPdq%uhXU&}k-7>yHH%BLCSI?_<04(7sy$X# zNee#J(^nI7A8+IilpE3HH@UXthmed*DZ{lW!qjCQ_E;}83p3iLTkwEIyE`E-<`KESLB7F?)vJ-_s*bh}%4$b=6H3!_w(fDA2s2?$MH3v#b2qw| z^dHvyg_~clE-l~|LL)P=a_hApMf>{V)?%`vDH&n|SQpdSxJEHbzI6t=yO`p0!-O_K zi}S9nxIg!C=$kg>TPCH)vw=A{tSi(y=^Dl3BDu@6F;gC^cQsf}MI|Tv(Ma(^^i30G zGlk8>0kmkP1qy-D$RO~D^hsHhEz0udm9vNfU8FHUT zkCG0WY>U`8`KgH^3EfmloSd_{u{IRqCoeVTq25?#U@WWxTUrt&5b<)2!-++ccP+lm${hz z`4*Sy3FUGR!V8Kx#Al~7PV2M{7mjj1(DrX>mHyc%_iC^(uPp0}|JB8|wQ)O1$h@VV z{gMjYAG?NX39K=Y3kD=2(&s-xEW}aHYw>3xTCKbuWLqAa-|E$8YhXf+jFro6c$5?~l~_;~3JoF6={sPs>Df#qO}9<2 z5U1q5JRo6S%fShVH#Ik@J%5Ufh5&z*SVW#K3{zFfU|Bs^{QMbBc}*cRnyg-ZT%W%hKMM)20XiiR40>+!QbZNpL%;;&AnrJ*X62B&Ox-H zA@5Qqujs>?nK<$EJVtM_qg|m3bVDk1{;D?5V?kWpqIsZ*^o&F#17h7WrIKOmcRI^s zd%UPwg<3ujMeY;3X_?$@M|PglZdvMG*%RCJSVO})Cu!5jh-Oz=HUolilKoj50(6(qRtfUA9&xO!DB?lZXUmkv#6{aTBj~Jw7dHVsr)KT(UvzK`A!KlOI_= zV2eh+L8gPHPN_XQ+#)CA)w@Q`2#=bI!V>y+b5{-J+EY?!)t(d@Mew`}jsmohE+>m`r!>LWzm#KA5;(7+ue0j1v16aX-n+;P*#QjZo15(61}AFIE%#}C{Yg#5u*o-!f2Q#{mtORs_*|lS zb6dD9>2H&TjIhX)Y$jCBd|kZsc=&R$?%)SyVB8-2Su~zg=NbDh)=Q{}Ir@*Bb=!Px z#}0D7qt(iC`F>b4+pvoV5MQ(BEq%YFJt!_dyCz?4J4(Ft(5;!VOp{D5vl{btr@Y}a zt36wMaMRBW_pRGq&b)iBvqmpvsNV#C?eokYC-bzQji_bm(YdO}!KK|y1DRha&5IJP zOAmS$`6Al$5J3ilwU|>Kd-?YV`J?F^h%8%!wW?3lSEAJ{ho3{(>m~#S9-o;V-8BL6 zKj82#UmPE^j8uL(E+4U3l{{{-bTPRW3p$2a+$W-5@Fg|PBI;vNm3^D_$H|AS(btn z)&84E9hq|PVT%&VXJ^WJ3pa+cZ>wSnXq$wpQBXc^l+U-4#t;04^z9y7S7F(xI84=8KM;-XA(WFju>y!{65yZ=4JnHfVCg;UApuSS=xap*}4)TkB&)jV+LeGHhNI zOu3@gi!*k~pO4P0eI{yvJGD%dHuD#eyG4o@L(>-Zx0ae5rBc@&?TDA7be7|EQI{Xq zkM^x$m}LD4%ER8h^^D7^g~6A6PH0`j=Z}vRuP2Xl94rtwOKBuOm?{=aomr2%y=dO! ztW$7yQ|q9=+Wkt8+o%K8A6;(0Fkg5GO3x$W$b&_tlZvE>_f7V@>k2fcGY2pYCMC^# zUr&YEEa1OZEe!OT*YA&$T!^gfqvx*$!-}j;eLM8IuP{ZaIlRf3*~?M3cUMIc+Y*O& zcec;N{Oj!!6Cakob1h}1DeYKSa&gov#|D>mKgwF_HZ@e?5d54jgo;U4V)s`%&{WEk zRvr5RqosrlNp4YG>1*+@cd1$9#>HU8$XAKgeKj6!I%Na+=rm8a}m>$I}Eo0a>b;%NH#pcBDFWNm^9>Z_uCnM6uW;IP`+s0^1%S|F9 zOy#NmQ+QJKoIOwV^vcEq=H>e1VngN#wH%~kyip)YV>H9TBK zi!IJq=MP>gs6<(qnwSH_py1Rml#fW`*e2^#qRmF>8d6M4>6h6JmRB6%uiME-Hx)LV z^@BghQ%!CUuVKvdQG&2vqDJxL7US9xW^|Rc>n#_rm4f|^Nv^07WRJypuAW+&cRX?4 z{&!Z@K+26c44W)Xw+qbR~6363w$WEKF zxQ5xLfdn2-{kF0DeycuSs(Lv}3m+ISTDS!6Np> zK#9S;h5CYKTJoY!_ZCGjvu}2^`sSZEM|oR^?C_~PK5c~iXUgORDzn_>o0=|Dior=X zN11lpV{(>wf4axu^FQ-mb@#?oXP!Frxaatk5r2M7MLCrGd~|EN`jd@ag^i6&8hJ9) z>z0cf1-oyxZ;?Q#FP1Tb+l`|9)E&(qUmEjOH;r}2t=aUxm#)7|wmrJ_aN2}if6f4h zgQvcnp6p0TQ8srO=6-~Yt+q}}52Va>=c#O94rL^s-~5W5Wg7Kzr%-|>o8`wzlM}_v zBWoM0r?#R~OYaX#p>HP}h+RnH)iHetYDa-Fh}}F`tlm}@_da7<hPQjIXc%;@nTV@)3b^*mr*^cknWzV#s3{BaKL@bUEjw*R0PW%1mew5j#;a6h>XW zq{1xAMe4Q=Zrq6|1TN=9k%z(+6PX^3yW{b@7OSrbeCIMKC6C|N?*CO%PhJ@WM&zi| z5!-7D?ffGwh=+GbQrbH_YfR4ZZuIYkyd<4EVMf53XP6Zo7Bp4pv6+E!=u%2QQ6!0q zQj1I>;7GrGJZ2D=@QSrjz%dfJ)8 z55v)U77CT<$CLE`N0!3>DJeU4gMnV3MD>>SJ`X9f!$Ea&LN>-P!q|a7&{*U;*XZYc zU<>%LK;=vQH#`fRQJ`_LvWnM00lih|9(#PizqB==@Ls6X1lG|u=bN6W?YZiE1bAI+ z8aE6LKkVkZU7+ZJdiC;OGMjF_ir{;3P4!Rn4|C|*NLUQvfNe4ZP_?^`Ii}N9d;iUC zrcu<;`gSZf>WN=m4a*&%{;zFkkUG!hNJ3ysTOER+X{l$Ul~RNQ^4gpscn~Qr7+o4F ze(#uuzU=CaE)wxf$l(__2hMTkjaDA`Bev8dFD9>{whUT7PLc2v5j)AQ*!a@2cU9@d?ODj8-i{!}f^ zIVkeT*cXpvHbS{GMj;o7;&hxj+RD5?FZhXUkY4+l_`F;HGG*t7BWC}t<0|#G^6(IP z)$RYZ6`ezBMBCoxy4&OalU~#586G)Y2D!8?s<)_d_^WnRf7exPm&YD~g-4LVE6>VI zps-I}ZivHfc&yJ$ClO~OQZPu!cle=CywP;SxrF6~R{;4P@nrtAlhOTdh%O3fwdSh2 zIQm-{8n$7T;iQ6U#j0q9zG($&9NxZL>?zQUaTTc8@F)~Bt$c=P6O;Wq z!Op5DTQiQS+r@v_aQ-SNHC!Y3G?{%j)_^Ax@9#+SvQ^Bcw2L2T}bCL$^~3Jcg+$238j>ZxKHDi)F55 zmjsr7THhc@=?z@FUpSGU-NJ6Jaxu$U9#gcuI{I$7sp)li3dB_y{2?d#`JR);jFnai{$8hn_mNG)sRPvgvBjS?^v;er;6-xC#a=e zEB&mL$4I&7p~(7s7WT9QQTkaAs7T1;0$fi?JqtChTxFjFUTzZ#AWs7M?pvsf|4?cS zRo`AP=ZM~?BN*n|3a%6+_gk|~^QH#(73hsAA3@$ex#%Oa-@&c|DCt0~Hd#|L2xbOf zr!7Jw+~(T)ZaGnF?VK|Dy}0sBKt1evCEEDwEs3?)o-2pB%1`7c*(VZoU5M`Tkn6@< zUNyZOSeryCoO=5fiai>NDRUo*#$Mi0t()EXQ|Os_^cd@@uVmxld{k?7e+WQPkL36()Ro zw-`F;$Vd9Mo?twOD;dZ)C>5E`yyYhxOtuv8%8#|&Jlu zbD@<1E~xY|9L*yzq)%mz&&{_l?mTziJ4Ra&mbB`cr;$8nb@FV^g&qlP|3wwFjgSWu z0v+&Cl~={}$gr0IB#uWVD^teL5wpLbv#wIaNdmu^amBg2xo(|T&}l2`K5jlrzo~w+ z{f`0_b-G^@R05tr)F{3dgZohFs=kBU{1h>}tDUQ+g&KU;a?I+Kj*>{?R!^kNvu9hg zmjh?!vz_klF!v8~XqGnd4U#z2QLXC0*`vD|J8T+}L>Ygi5pf^p<(wK_Fl&qW8`nBvG*966BqnyuRJt+u`wwOI>oG7aa5Sy1iE!K zVBC$b29IWKvlaSFM;l<&6vq=;ri`Qmgss=bX$tcOw9G%sgkfLNc_14=eS1~rCKeAw zkz{979JL7g7l^heG`ruB`DZ4A-E)I%$OSbQ;x#bCfR!2z`~bhlMa4(da9AjWvzWW}NzJ`fU za=F8@@Ad_d<-Sa*&C0u)8HS@b@EoUNwllZXhR8i!S;v|90Q4;Qx%Ie5XmJ#7r+}l= zC8u+C$TH7^qAiVT8&oe^> zF*@RF)#J_{(}%0F`Guvw*(viMrh7}^wwe3?Lx%a!RW`TMl^aVUTwP-6m)FwJRnNF& H`RsoH%K*So literal 0 HcmV?d00001 From 9b67654810eaa4e61d2a8d09d3e669202e4c2635 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Sun, 31 May 2026 21:27:18 -0500 Subject: [PATCH 2/3] Implement ESDC1A exciter model --- CHANGELOG.md | 1 + .../Model/PhasorDynamics/ComponentLibrary.hpp | 1 + .../PhasorDynamics/Exciter/CMakeLists.txt | 1 + .../Exciter/ESDC1A/CMakeLists.txt | 49 ++ .../PhasorDynamics/Exciter/ESDC1A/Esdc1a.cpp | 27 + .../PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp | 152 +++++ .../Exciter/ESDC1A/Esdc1aData.hpp | 76 +++ .../ESDC1A/Esdc1aDependencyTracking.cpp | 27 + .../Exciter/ESDC1A/Esdc1aEnzyme.cpp | 93 +++ .../Exciter/ESDC1A/Esdc1aImpl.hpp | 480 ++++++++++++++ .../PhasorDynamics/Exciter/ESDC1A/README.md | 45 +- .../Model/PhasorDynamics/Exciter/README.md | 2 +- GridKit/Model/PhasorDynamics/INPUT_FORMAT.md | 1 + GridKit/Model/PhasorDynamics/SystemModel.hpp | 37 ++ .../Model/PhasorDynamics/SystemModelData.hpp | 3 + .../SystemModelDataJSONParser.hpp | 6 + tests/UnitTests/PhasorDynamics/CMakeLists.txt | 8 + .../PhasorDynamics/ExciterEsdc1aTests.hpp | 585 ++++++++++++++++++ .../PhasorDynamics/runExciterEsdc1aTests.cpp | 18 + 19 files changed, 1591 insertions(+), 21 deletions(-) create mode 100644 GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt create mode 100644 GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.cpp create mode 100644 GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp create mode 100644 GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp create mode 100644 GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp create mode 100644 GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp create mode 100644 GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp create mode 100644 tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp create mode 100644 tests/UnitTests/PhasorDynamics/runExciterEsdc1aTests.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index 0c31a6e41..0fdc82a38 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -51,6 +51,7 @@ - Added component model developer checklist to a README file. - Added `IEEEST` Stabilizer Model - Added `SEXS-PTI` Exciter Model +- Added `ESDC1A` Exciter Model - Added `GENSAL` Machine Model - Added 200 Bus Synthetic Illinois Case - Added node objects to `PowerElectronics` module & updated all examples to make use of them. diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index e48a8a9c5..84cd57961 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Exciter/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Exciter/CMakeLists.txt index 1120e20c2..7e5e4b9e3 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Exciter/CMakeLists.txt @@ -3,5 +3,6 @@ # - Luke Lowery # ]] +add_subdirectory(ESDC1A) add_subdirectory(IEEET1) add_subdirectory(SEXS-PTI) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt new file mode 100644 index 000000000..6886a26ec --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt @@ -0,0 +1,49 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers + Esdc1a.hpp + Esdc1aData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library(phasor_dynamics_exciter_esdc1a + SOURCES + Esdc1aEnzyme.cpp + HEADERS + ${_install_headers} + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core + PUBLIC GridKit::phasor_dynamics_signal + PRIVATE ClangEnzymeFlags + COMPILE_OPTIONS + PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) +else() + gridkit_add_library(phasor_dynamics_exciter_esdc1a + SOURCES + Esdc1a.cpp + HEADERS + ${_install_headers} + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + GridKit::phasor_dynamics_core + GridKit::phasor_dynamics_signal) +endif() + +gridkit_add_library(phasor_dynamics_exciter_esdc1a_dependency_tracking + SOURCES + Esdc1aDependencyTracking.cpp + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + GridKit::phasor_dynamics_core + GridKit::phasor_dynamics_signal_dependency_tracking) + +target_link_libraries(phasor_dynamics_components + INTERFACE GridKit::phasor_dynamics_exciter_esdc1a) +target_link_libraries(phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_exciter_esdc1a_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.cpp new file mode 100644 index 000000000..6c5c0f89c --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.cpp @@ -0,0 +1,27 @@ +/** + * @file Esdc1a.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the ESDC1A exciter model. + */ + +#include "Esdc1aImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + template + int Esdc1a::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Esdc1a..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + return 0; + } + + template class Esdc1a; + template class Esdc1a; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp new file mode 100644 index 000000000..d401eac00 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp @@ -0,0 +1,152 @@ +/** + * @file Esdc1a.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the ESDC1A exciter model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + template + class SignalNode; + + namespace Exciter + { + /// Internal variables of an `Esdc1a`. + enum class Esdc1aInternalVariables : size_t + { + EFDP, ///< Field-voltage state before optional speed multiplier + VC, ///< Sensed compensated voltage + VR, ///< Voltage-regulator output + VF, ///< Stabilizing feedback output + XLL, ///< Lead-lag state + EV, ///< Voltage-regulator input error + VLL, ///< Lead-lag block output + VHV, ///< High-value gate output + SE, ///< Saturation coefficient + VFE, ///< Exciter feedback signal + EFD, ///< Field-voltage output + MAXIMUM, + }; + + /// External variables of an `Esdc1a`. + enum class Esdc1aExternalVariables : size_t + { + OMEGA, ///< Machine speed deviation + VS, ///< Stabilizer input signal + VUEL, ///< Under-excitation limiter input + MAXIMUM, + }; + + template + class Esdc1a : public Component + { + using Component::alpha_; + using Component::f_; + using Component::gridkit_component_id_; + using Component::J_; + using Component::J_cols_buffer_; + using Component::J_rows_buffer_; + using Component::J_vals_buffer_; + using Component::residual_indices_; + using Component::size_; + using Component::tag_; + using Component::variable_indices_; + using Component::wb_; + using Component::y_; + using Component::yp_; + + public: + using RealT = typename Component::RealT; + using bus_type = BusBase; + using signal_type = SignalNode; + using model_data_type = Esdc1aData; + using MonitorT = Model::VariableMonitor; + + Esdc1a(bus_type* bus); + Esdc1a(bus_type* bus, const model_data_type& data); + ~Esdc1a(); + + int setGridKitComponentID(IdxT) override final; + int allocate() override final; + int verify() const override final; + int initialize() override final; + int tagDifferentiable() override final; + int evaluateResidual() override final; + int evaluateJacobian() override final; + + auto getSignals() + -> ComponentSignals& + { + return signals_; + } + + const Model::VariableMonitorBase* getMonitor() const override; + + __attribute__((always_inline)) inline int evaluateInternalResidual( + ScalarT*, ScalarT*, ScalarT*, ScalarT*, ScalarT*); + + private: + void initModelParams(const model_data_type& data); + void setDerivedParams(); + void initializeMonitor(); + + bus_type* bus_{nullptr}; + + RealT Tr_{0.0}; + RealT Ka_{40.0}; + RealT Ta_{0.1}; + RealT Tb_{0.0}; + RealT Tc_{0.0}; + RealT Vrmax_{1.0}; + RealT Vrmin_{-1.0}; + RealT Ke_{0.1}; + RealT Te_{0.5}; + RealT Kf_{0.05}; + RealT Tf1_{0.7}; + RealT spdmlt_{0.0}; + RealT E1_{2.8}; + RealT Se1_{0.08}; + RealT E2_{3.7}; + RealT Se2_{0.33}; + IdxT UEL_{0}; + RealT exclim_{1.0}; + + RealT sUEL_{0}; + RealT sUELoff_{1}; + RealT slim_{0}; + RealT slim_off_{1}; + RealT use_lead_lag_{0}; + RealT bypass_lead_lag_{1}; + RealT SA_{0}; + RealT SB_{0}; + + ScalarT vref_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp new file mode 100644 index 000000000..0c6a23797 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp @@ -0,0 +1,76 @@ +/** + * @file Esdc1aData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the ESDC1A exciter model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + /// Parameter keys for the ESDC1A exciter model. + enum class Esdc1aParameters + { + Tr, ///< Transducer time constant + Ka, ///< Voltage-regulator gain + Ta, ///< Voltage-regulator time constant + Tb, ///< Lead-lag denominator time constant + Tc, ///< Lead-lag numerator time constant + Vrmax, ///< Maximum voltage-regulator output + Vrmin, ///< Minimum voltage-regulator output + Ke, ///< Exciter field-resistance line-slope margin + Te, ///< Exciter field time constant + Kf, ///< Stabilizing feedback gain + Tf1, ///< Feedback lead time constant + Spdmlt, ///< Speed multiplier flag + E1, ///< First saturation voltage point + Se1, ///< Saturation value at E1 + E2, ///< Second saturation voltage point + Se2, ///< Saturation value at E2 + UEL, ///< UEL input-location selector + exclim ///< Exciter feedback lower-limit flag + }; + + /// Ports for the ESDC1A exciter model. + enum class Esdc1aPorts + { + bus, ///< Unique ID of the terminal bus + speed, ///< Unique ID of the generator speed-deviation signal + efd, ///< Unique ID of the output EFD signal + vs, ///< Unique ID of the optional stabilizer input signal + vuel ///< Unique ID of the optional UEL input signal + }; + + /// Variables available through the monitor interface. + enum class Esdc1aMonitorableVariables + { + efd, ///< Field-voltage output + vc, ///< Sensed compensated voltage + vr, ///< Voltage-regulator output + vf, ///< Stabilizing feedback state + se, ///< Saturation coefficient + vfe ///< Exciter feedback signal + }; + + template + struct Esdc1aData : public ComponentData + { + Esdc1aData() = default; + + using Parameters = Esdc1aParameters; + using Ports = Esdc1aPorts; + using MonitorableVariables = Esdc1aMonitorableVariables; + }; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp new file mode 100644 index 000000000..c7c25aec6 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file Esdc1aDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the ESDC1A exciter model. + */ + +#include "Esdc1aImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + template + int Esdc1a::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Esdc1a..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + return 0; + } + + template class Esdc1a; + template class Esdc1a; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp new file mode 100644 index 000000000..96a7837bd --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp @@ -0,0 +1,93 @@ +/** + * @file Esdc1aEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the ESDC1A exciter model. + */ + +#include + +#include "Esdc1aImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + template + int Esdc1a::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Esdc1a..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); + if (J_rows_buffer_ == nullptr) + { + J_rows_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; + J_cols_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; + J_vals_buffer_ = new RealT[static_cast(size_) * static_cast(size_)]; + } + + using ModelT = GridKit::PhasorDynamics::Exciter::Esdc1a; + using Fn = GridKit::Enzyme::Sparse::MemberFunctions; + + GridKit::Enzyme::Sparse::DfDy::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + GridKit::Enzyme::Sparse::DfDwb::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + GridKit::Enzyme::Sparse::DfDws::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + return 0; + } + + template class Esdc1a; + template class Esdc1a; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp new file mode 100644 index 000000000..87d648eb6 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -0,0 +1,480 @@ +/** + * @file Esdc1aImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the ESDC1A exciter model. + */ + +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + using Log = ::GridKit::Utilities::Logger; + + template + Esdc1a::Esdc1a(bus_type* bus) + : bus_(bus) + { + setDerivedParams(); + size_ = static_cast(Esdc1aInternalVariables::MAXIMUM); + } + + template + Esdc1a::Esdc1a(bus_type* bus, const model_data_type& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initModelParams(data); + setDerivedParams(); + initializeMonitor(); + size_ = static_cast(Esdc1aInternalVariables::MAXIMUM); + } + + template + Esdc1a::~Esdc1a() + { + } + + template + void Esdc1a::initModelParams(const model_data_type& data) + { + using Params = typename model_data_type::Parameters; + + auto loadReal = [&](auto param, RealT& member) + { + if (data.parameters.contains(param)) + { + member = std::get(data.parameters.at(param)); + } + }; + + auto loadIndex = [&](auto param, IdxT& member) + { + if (data.parameters.contains(param)) + { + member = std::get(data.parameters.at(param)); + } + }; + + loadReal(Params::Tr, Tr_); + loadReal(Params::Ka, Ka_); + loadReal(Params::Ta, Ta_); + loadReal(Params::Tb, Tb_); + loadReal(Params::Tc, Tc_); + loadReal(Params::Vrmax, Vrmax_); + loadReal(Params::Vrmin, Vrmin_); + loadReal(Params::Ke, Ke_); + loadReal(Params::Te, Te_); + loadReal(Params::Kf, Kf_); + loadReal(Params::Tf1, Tf1_); + loadReal(Params::Spdmlt, spdmlt_); + loadReal(Params::E1, E1_); + loadReal(Params::Se1, Se1_); + loadReal(Params::E2, E2_); + loadReal(Params::Se2, Se2_); + loadIndex(Params::UEL, UEL_); + loadReal(Params::exclim, exclim_); + } + + template + void Esdc1a::setDerivedParams() + { + sUEL_ = ZERO; + if (UEL_ >= static_cast(2)) + { + sUEL_ = ONE; + } + sUELoff_ = ONE - sUEL_; + slim_ = exclim_; + slim_off_ = ONE - slim_; + + use_lead_lag_ = ZERO; + if (Tb_ > ZERO) + { + use_lead_lag_ = ONE; + } + bypass_lead_lag_ = ONE - use_lead_lag_; + + if (Se1_ == ZERO && Se2_ == ZERO) + { + SA_ = ZERO; + SB_ = ZERO; + return; + } + if (E1_ <= ZERO || E2_ <= ZERO || E1_ == E2_ + || Se1_ <= ZERO || Se2_ <= ZERO || Se1_ == Se2_) + { + SA_ = ZERO; + SB_ = ZERO; + return; + } + + const RealT C = std::sqrt(Se2_ / Se1_); + SA_ = (C * E1_ - E2_) / (C - ONE); + SB_ = Se1_ / ((E1_ - SA_) * (E1_ - SA_)); + } + + template + const Model::VariableMonitorBase* Esdc1a::getMonitor() const + { + return monitor_.get(); + } + + template + void Esdc1a::initializeMonitor() + { + using Variable = typename model_data_type::MonitorableVariables; + auto index = [](Esdc1aInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::efd, [this, index] + { return y_[index(Esdc1aInternalVariables::EFD)]; }); + monitor_->set(Variable::vc, [this, index] + { return y_[index(Esdc1aInternalVariables::VC)]; }); + monitor_->set(Variable::vr, [this, index] + { return y_[index(Esdc1aInternalVariables::VR)]; }); + monitor_->set(Variable::vf, [this, index] + { return y_[index(Esdc1aInternalVariables::VF)]; }); + monitor_->set(Variable::se, [this, index] + { return y_[index(Esdc1aInternalVariables::SE)]; }); + monitor_->set(Variable::vfe, [this, index] + { return y_[index(Esdc1aInternalVariables::VFE)]; }); + } + + template + int Esdc1a::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Esdc1a::allocate() + { + size_ = static_cast(Esdc1aInternalVariables::MAXIMUM); + auto size = static_cast(size_); + + f_.assign(size, ScalarT{0}); + y_.assign(size, ScalarT{0}); + yp_.assign(size, ScalarT{0}); + tag_.assign(size, false); + variable_indices_.resize(size); + residual_indices_.resize(size); + + wb_.assign(2, ScalarT{0}); + + auto signal_size = static_cast(Esdc1aExternalVariables::MAXIMUM); + ws_.assign(signal_size, ScalarT{0}); + ws_indices_.assign(signal_size, INVALID_INDEX); + + for (IdxT j = 0; j < size_; ++j) + { + this->setVariableIndex(j, j); + this->setResidualIndex(j, j); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(Esdc1aInternalVariables::EFD)], + &(this->getVariableIndex(static_cast(Esdc1aInternalVariables::EFD)))); + } + + return 0; + } + + template + int Esdc1a::verify() const + { + int ret = 0; + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Esdc1a: " << message << '\n'; + ret += 1; + } + }; + + if (bus_ == nullptr) + { + Log::error() << "Esdc1a: bus pointer is null\n"; + ret += 1; + } + + check(Ka_ > ZERO, "Ka must be positive"); + check(Tr_ >= ZERO, "Tr must be non-negative"); + check(Ta_ > ZERO, "Ta must be positive"); + check(Tb_ >= ZERO, "Tb must be non-negative"); + check(Tc_ >= ZERO, "Tc must be non-negative"); + check(Tb_ > ZERO || Tc_ == ZERO, "Tc must be zero when Tb is zero"); + check(Te_ > ZERO, "Te must be positive"); + check(Tf1_ >= ZERO, "Tf1 must be non-negative"); + check(Vrmin_ <= Vrmax_, "Vrmin must be less than or equal to Vrmax"); + check(spdmlt_ == ZERO || spdmlt_ == ONE, "Spdmlt must be 0 or 1"); + check(exclim_ == ZERO || exclim_ == ONE, "exclim must be 0 or 1"); + check(static_cast(UEL_) >= ZERO && UEL_ <= static_cast(3), + "UEL must be 0, 1, 2, or 3"); + + if (Se1_ == ZERO && Se2_ == ZERO) + { + // Saturation disabled. + } + else + { + check(E1_ > ZERO, "E1 must be positive when saturation is enabled"); + check(E2_ > ZERO, "E2 must be positive when saturation is enabled"); + check(E1_ != E2_, "E1 and E2 must differ when saturation is enabled"); + check(Se1_ > ZERO, "Se1 must be positive when saturation is enabled"); + check(Se2_ > ZERO, "Se2 must be positive when saturation is enabled"); + check(Se1_ != Se2_, "Se1 and Se2 must differ when saturation is enabled"); + } + + if (!signals_.template isAssigned()) + { + Log::error() << "Esdc1a: required EFD signal is not assigned\n"; + ret += 1; + } + + if (spdmlt_ == ONE + && !signals_.template isAttached()) + { + Log::error() << "Esdc1a: speed signal is required when Spdmlt is enabled\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Esdc1a: speed signal attached with no linked source\n"; + ret += 1; + } + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Esdc1a: VS signal attached with no linked source\n"; + ret += 1; + } + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Esdc1a: VUEL signal attached with no linked source\n"; + ret += 1; + } + + return ret; + } + + template + int Esdc1a::initialize() + { + const auto EFDP = static_cast(Esdc1aInternalVariables::EFDP); + const auto VC = static_cast(Esdc1aInternalVariables::VC); + const auto VR = static_cast(Esdc1aInternalVariables::VR); + const auto VF = static_cast(Esdc1aInternalVariables::VF); + const auto XLL = static_cast(Esdc1aInternalVariables::XLL); + const auto EV = static_cast(Esdc1aInternalVariables::EV); + const auto VLL = static_cast(Esdc1aInternalVariables::VLL); + const auto VHV = static_cast(Esdc1aInternalVariables::VHV); + const auto SE = static_cast(Esdc1aInternalVariables::SE); + const auto VFE = static_cast(Esdc1aInternalVariables::VFE); + const auto EFD = static_cast(Esdc1aInternalVariables::EFD); + + ScalarT omega0{ZERO}; + if (signals_.template isAttached()) + { + omega0 = signals_.template readExternalVariable(); + } + + ScalarT vs0{ZERO}; + if (signals_.template isAttached()) + { + vs0 = signals_.template readExternalVariable(); + } + + ScalarT vuel0{ZERO}; + if (signals_.template isAttached()) + { + vuel0 = signals_.template readExternalVariable(); + } + + const ScalarT denom = ONE + spdmlt_ * omega0; + if (denom == ZERO) + { + Log::error() << "Esdc1a: speed multiplier denominator is zero at initialization\n"; + return 1; + } + + const ScalarT Ec0 = std::sqrt(bus_->Vr() * bus_->Vr() + bus_->Vi() * bus_->Vi()); + + const ScalarT efd0 = y_[EFD]; + const ScalarT efdp0 = efd0 / denom; + const ScalarT se0 = SB_ * Math::qramp(efdp0 - SA_); + const ScalarT vfe0 = slim_off_ * (Ke_ + se0) * efdp0 + + slim_ * Math::ramp((Ke_ + se0) * efdp0); + const ScalarT vr0 = vfe0; + const ScalarT vhv0 = vr0 / Ka_; + const ScalarT vc0 = Ec0; + const ScalarT vf0 = ScalarT{ZERO}; + const ScalarT xll0 = vhv0; + const ScalarT vll0 = vhv0; + const ScalarT ev0 = vhv0; + + if (vr0 < Vrmin_ || vr0 > Vrmax_) + { + Log::error() << "Esdc1a: initialized VR is outside limits\n"; + return 1; + } + if (sUEL_ == ZERO && vhv0 < vuel0) + { + Log::error() << "Esdc1a: high-value gate active at initialization\n"; + return 1; + } + + vref_ = ev0 + vc0 + vf0 - vs0 - sUEL_ * vuel0; + + y_[EFDP] = efdp0; + y_[VC] = vc0; + y_[VR] = vr0; + y_[VF] = vf0; + y_[XLL] = xll0; + y_[EV] = ev0; + y_[VLL] = vll0; + y_[VHV] = vhv0; + y_[SE] = se0; + y_[VFE] = vfe0; + y_[EFD] = efd0; + + std::fill(yp_.begin(), yp_.end(), ZERO); + return 0; + } + + template + int Esdc1a::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(Esdc1aInternalVariables::EFDP)] = true; + tag_[static_cast(Esdc1aInternalVariables::VC)] = (Tr_ > ZERO); + tag_[static_cast(Esdc1aInternalVariables::VR)] = true; + tag_[static_cast(Esdc1aInternalVariables::VF)] = (Tf1_ > ZERO); + tag_[static_cast(Esdc1aInternalVariables::XLL)] = (Tb_ > ZERO); + return 0; + } + + template + __attribute__((always_inline)) inline int Esdc1a::evaluateInternalResidual( + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + ScalarT* f) + { + const auto EFDP = static_cast(Esdc1aInternalVariables::EFDP); + const auto VC = static_cast(Esdc1aInternalVariables::VC); + const auto VR = static_cast(Esdc1aInternalVariables::VR); + const auto VF = static_cast(Esdc1aInternalVariables::VF); + const auto XLL = static_cast(Esdc1aInternalVariables::XLL); + const auto EV = static_cast(Esdc1aInternalVariables::EV); + const auto VLL = static_cast(Esdc1aInternalVariables::VLL); + const auto VHV = static_cast(Esdc1aInternalVariables::VHV); + const auto SE = static_cast(Esdc1aInternalVariables::SE); + const auto VFE = static_cast(Esdc1aInternalVariables::VFE); + const auto EFD = static_cast(Esdc1aInternalVariables::EFD); + + const auto OMEGA = static_cast(Esdc1aExternalVariables::OMEGA); + const auto VS = static_cast(Esdc1aExternalVariables::VS); + const auto VUEL = static_cast(Esdc1aExternalVariables::VUEL); + + const ScalarT efdp = y[EFDP]; + const ScalarT vc = y[VC]; + const ScalarT vr = y[VR]; + const ScalarT vf = y[VF]; + const ScalarT xll = y[XLL]; + const ScalarT ev = y[EV]; + const ScalarT vll = y[VLL]; + const ScalarT vhv = y[VHV]; + const ScalarT se = y[SE]; + const ScalarT vfe = y[VFE]; + const ScalarT efd = y[EFD]; + + const ScalarT omega = ws[OMEGA]; + const ScalarT vs = ws[VS]; + const ScalarT vuel = ws[VUEL]; + + const ScalarT Ec = std::sqrt(wb[0] * wb[0] + wb[1] * wb[1]); + + f[EFDP] = -Te_ * yp[EFDP] + vr - vfe; + f[VC] = -Tr_ * yp[VC] - vc + Ec; + f[VR] = -Ta_ * yp[VR] + Math::antiwindup(vr, -vr + Ka_ * vhv, Vrmin_, Vrmax_); + f[VF] = -Te_ * Tf1_ * yp[VF] - Te_ * vf + Kf_ * (vr - vfe); + f[XLL] = -Tb_ * yp[XLL] - xll + ev; + + // Keep the EV equation grouped this way to avoid an Enzyme sparse-lowering compile failure. + const ScalarT ref_err = vref_ - ev; + const ScalarT sig_err = vs - vc; + const ScalarT lim_err = sUEL_ * vuel - vf; + f[EV] = ref_err + sig_err + lim_err; + f[VLL] = use_lead_lag_ * (-Tb_ * (vll - xll) + Tc_ * (ev - xll)) + + bypass_lead_lag_ * (-vll + ev); + f[VHV] = -vhv + sUEL_ * vll + sUELoff_ * Math::max(vll, vuel); + f[SE] = -se + SB_ * Math::qramp(efdp - SA_); + f[VFE] = -vfe + slim_off_ * (Ke_ + se) * efdp + slim_ * Math::ramp((Ke_ + se) * efdp); + f[EFD] = -efd + (ONE + spdmlt_ * omega) * efdp; + + return 0; + } + + template + int Esdc1a::evaluateResidual() + { + const auto OMEGA = static_cast(Esdc1aExternalVariables::OMEGA); + const auto VS = static_cast(Esdc1aExternalVariables::VS); + const auto VUEL = static_cast(Esdc1aExternalVariables::VUEL); + + std::fill(ws_.begin(), ws_.end(), ScalarT{ZERO}); + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + if (signals_.template isAttached()) + { + ws_[OMEGA] = signals_.template readExternalVariable(); + ws_indices_[OMEGA] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[VS] = signals_.template readExternalVariable(); + ws_indices_[VS] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[VUEL] = signals_.template readExternalVariable(); + ws_indices_[VUEL] = + signals_.template readExternalVariableIndex(); + } + + wb_[0] = bus_->Vr(); + wb_[1] = bus_->Vi(); + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + return 0; + } + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md index e70e9c74e..1161ff2a7 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md @@ -19,24 +19,24 @@ Notes: Symbol | Units | Description | Typical Value | Note ------------------------------------|----------|---------------------------------------------------------|---------------|------ -$T_R$ | [sec] | Transducer time constant | TBD | Block name: `Tr`; if zero, $V_C$ is algebraic -$K_A$ | [p.u.] | Voltage-regulator gain | TBD | Block name: `Ka` -$T_A$ | [sec] | Voltage-regulator time constant | TBD | Block name: `Ta` -$T_B$ | [sec] | Lag time constant for the voltage-regulator input lead-lag block | TBD | Block name: `Tb` -$T_C$ | [sec] | Lead time constant for the voltage-regulator input lead-lag block | TBD | Block name: `Tc` -$V_R^{\max}$ | [p.u.] | Maximum voltage-regulator output | TBD | Block name: `Vrmax` -$V_R^{\min}$ | [p.u.] | Minimum voltage-regulator output | TBD | Block name: `Vrmin` -$K_E$ | [p.u.] | Exciter field-resistance line-slope margin | TBD | Block name: `Ke` -$T_E$ | [sec] | Exciter time constant | TBD | Block name: `Te` -$K_F$ | [p.u.] | Stabilizing feedback gain | TBD | Block name: `Kf` -$T_{F1}$ | [sec] | Feedback lead time constant | TBD | Block name: `Tf1` -$s_{\mathrm{spd}}$ | [binary] | Speed multiplier flag | TBD | Block name: `Spdmlt`; nonzero source values are enabled -$E_1$ | [p.u.] | First saturation voltage point | TBD | Block name: `E1` -$S_E(E_1)$ | [p.u.] | Saturation value at $E_1$ | TBD | Block name: `Se1` -$E_2$ | [p.u.] | Second saturation voltage point | TBD | Block name: `E2` -$S_E(E_2)$ | [p.u.] | Saturation value at $E_2$ | TBD | Block name: `Se2` -$I_{\mathrm{UEL}}$ | [integer] | Under-excitation limiter input-location selector | TBD | Block name: `UEL`; 0/1 = HV gate input, 2/3 = input-error summing junction -$s_{\mathrm{lim}}$ | [binary] | Exciter feedback lower-limit flag | TBD | Block name: `exclim`; nonzero enables zero lower limit on $V_{\mathrm{FE}}$ +$T_R$ | [sec] | Transducer time constant | 0 | Block name: `Tr`; if zero, $V_C$ is algebraic +$K_A$ | [p.u.] | Voltage-regulator gain | 40 | Block name: `Ka` +$T_A$ | [sec] | Voltage-regulator time constant | 0.1 | Block name: `Ta` +$T_B$ | [sec] | Lag time constant for the voltage-regulator input lead-lag block | 0 | Block name: `Tb`; if $T_B=T_C=0$, GridKit treats the lead-lag block as bypassed +$T_C$ | [sec] | Lead time constant for the voltage-regulator input lead-lag block | 0 | Block name: `Tc`; must be zero when $T_B=0$ +$V_R^{\max}$ | [p.u.] | Maximum voltage-regulator output | 1 | Block name: `Vrmax` +$V_R^{\min}$ | [p.u.] | Minimum voltage-regulator output | -1 | Block name: `Vrmin` +$K_E$ | [p.u.] | Exciter field-resistance line-slope margin | 0.1 | Block name: `Ke` +$T_E$ | [sec] | Exciter time constant | 0.5 | Block name: `Te` +$K_F$ | [p.u.] | Stabilizing feedback gain | 0.05 | Block name: `Kf` +$T_{F1}$ | [sec] | Feedback lead time constant | 0.7 | Block name: `Tf1` +$s_{\mathrm{spd}}$ | [binary] | Speed multiplier flag | 0 | Block name: `Spdmlt`; 1 enables the speed multiplier +$E_1$ | [p.u.] | First saturation voltage point | 2.8 | Block name: `E1` +$S_E(E_1)$ | [p.u.] | Saturation value at $E_1$ | 0.08 | Block name: `Se1` +$E_2$ | [p.u.] | Second saturation voltage point | 3.7 | Block name: `E2` +$S_E(E_2)$ | [p.u.] | Saturation value at $E_2$ | 0.33 | Block name: `Se2` +$I_{\mathrm{UEL}}$ | [integer] | Under-excitation limiter input-location selector | 0 | Block name: `UEL`; 0/1 = HV gate input, 2/3 = input-error summing junction +$s_{\mathrm{lim}}$ | [binary] | Exciter feedback lower-limit flag | 1 | Block name: `exclim`; 1 enables the zero lower limit on $V_{\mathrm{FE}}$ ### Parameter Validation @@ -47,7 +47,8 @@ The required checks are: ```math \begin{aligned} &K_A > 0 \\ - &T_R \ge 0,\quad T_A > 0,\quad T_B > 0,\quad T_C \ge 0,\quad T_E > 0,\quad T_{F1} \ge 0 \\ + &T_R \ge 0,\quad T_A > 0,\quad T_B \ge 0,\quad T_C \ge 0,\quad T_E > 0,\quad T_{F1} \ge 0 \\ + &T_B > 0\quad\text{or}\quad(T_B = 0\ \text{and}\ T_C = 0) \\ &V_R^{\min} \le V_R^{\max} \\ &s_{\mathrm{spd}}, s_{\mathrm{lim}} \in \{0,1\} \\ &I_{\mathrm{UEL}} \in \{0,1,2,3\} @@ -87,6 +88,8 @@ The UEL routing flag and off-mode flag complements are: \end{aligned} ``` +The PowerWorld default $T_B=T_C=0$ is treated as a lead-lag bypass. + The saturation curve is fitted from the two supplied saturation points. If both saturation factors are zero, use $S_A=0$ and $S_B=0$. Otherwise: ```math @@ -167,7 +170,7 @@ The algebraic targets use CommonMath helper notation where applicable: ```math \begin{aligned} 0 &= -e_V + V_{\mathrm{ref}} + V_S + s_{\mathrm{UEL}}V_{\mathrm{UEL}} - V_C - V_F \\ - 0 &= -T_B V_{\mathrm{LL}} + T_C e_V + (T_B - T_C)x_{\mathrm{LL}} \\ + 0 &= -T_B(V_{\mathrm{LL}} - x_{\mathrm{LL}}) + T_C(e_V - x_{\mathrm{LL}}) \\ 0 &= -V_{\mathrm{HV}} + s_{\mathrm{UEL}}V_{\mathrm{LL}} + s_{\mathrm{UEL}}^{\mathrm{off}}\text{max}(V_{\mathrm{LL}}, V_{\mathrm{UEL}}) \\ @@ -181,6 +184,8 @@ The algebraic targets use CommonMath helper notation where applicable: CommonMath defines the helper targets and smooth approximations for [max](../../../../CommonMath.md#derived-functions) and the primitives [ramp and quadratic ramp](../../../../CommonMath.md#primitives) $\rho$ and $q$. +When the PowerWorld default $T_B=T_C=0$ is used, GridKit bypasses the lead-lag block so $V_{\mathrm{LL}}=e_V$. + ## Initialization The machine initializes $E_{\mathrm{fd}}$ first. For a standard unsaturated start, ESDC1A reads that value along with attached $\omega$, $E_C$, $V_S$, and $V_{\mathrm{UEL}}$, sets all internal derivatives to zero, and evaluates: diff --git a/GridKit/Model/PhasorDynamics/Exciter/README.md b/GridKit/Model/PhasorDynamics/Exciter/README.md index f204ca63f..e0ccf716c 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/README.md @@ -1,7 +1,7 @@ # **Exciter Models** > [!NOTE] -> IEEET1 and SEXS-PTI exciters are currently implemented. +> EXDC1 is not currently implemented. ## Introduction diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 32df36cdc..c7499f6f0 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -146,6 +146,7 @@ are specified: `GenClassical` | the classical machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Xdp`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega` `Tgov1 ` | the TGOV1 governor model | `pmech`, `speed` | `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` `Ieeet1` | the IEEET1 exciter model | `bus`, `speed`, `efd`, `vs`\* | `Tr`, `Ka`, `Ta`, `Ke`, `Te`, `Kf`, `Tf`, `Vrmin`, `Vrmax`, `E1`, `E2`, `Se1`, `Se2`, `Ispdlim` | `efd`, `ksat` + `Esdc1a` | the ESDC1A exciter model | `bus`, `efd`, `speed`\*, `vs`\*, `vuel`\* | `Tr`, `Ka`, `Ta`, `Tb`, `Tc`, `Vrmax`, `Vrmin`, `Ke`, `Te`, `Kf`, `Tf1`, `Spdmlt`, `E1`, `Se1`, `E2`, `Se2`, `UEL`, `exclim` | `efd`, `vc`, `vr`, `vf`, `se`, `vfe` `SexsPti` | the SEXS-PTI simplified exciter model | `bus`, `efd`, `vs`\* | `Ta`, `Tb`, `Te`, `K`, `Efdmax`, `Efdmin` | `efd` `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` `BusFault` | simple impedance-based fault at a bus | `bus`, `status`\* | `state0`, `R`, `X` | `state`, `ir`, `ii` diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index d7b4cc9d7..1d226cfb5 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -328,6 +328,43 @@ namespace GridKit addComponent(exciter); } + for (const auto& excitedata : data.esdc1a) + { + IdxT bus_index = 0; + if (excitedata.ports.contains(Esdc1aData::Ports::bus)) + { + bus_index = excitedata.ports.at(Esdc1aData::Ports::bus); + } + + auto* exciter = new Esdc1a(getBus(bus_index), excitedata); + + if (excitedata.ports.contains(Esdc1aData::Ports::speed)) + { + IdxT speed = excitedata.ports.at(Esdc1aData::Ports::speed); + exciter->getSignals().template attachSignalNode(getSignal(speed)); + } + + if (excitedata.ports.contains(Esdc1aData::Ports::efd)) + { + IdxT efd = excitedata.ports.at(Esdc1aData::Ports::efd); + exciter->getSignals().template assignSignalNode(getSignal(efd)); + } + + if (excitedata.ports.contains(Esdc1aData::Ports::vs)) + { + IdxT vs = excitedata.ports.at(Esdc1aData::Ports::vs); + exciter->getSignals().template attachSignalNode(getSignal(vs)); + } + + if (excitedata.ports.contains(Esdc1aData::Ports::vuel)) + { + IdxT vuel = excitedata.ports.at(Esdc1aData::Ports::vuel); + exciter->getSignals().template attachSignalNode(getSignal(vuel)); + } + + addComponent(exciter); + } + for (const auto& excitedata : data.sexspti) { IdxT bus_index = 0; diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index d4deef66e..67b765c20 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -41,6 +42,7 @@ namespace GridKit using BusToSignalAdapterDataT = BusToSignalAdapterData; using BusFaultDataT = BusFaultData; using Tgov1DataT = Governor::Tgov1Data; + using Esdc1aDataT = Exciter::Esdc1aData; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; using IeeestDataT = Stabilizer::IeeestData; @@ -99,6 +101,7 @@ namespace GridKit std::vector load; ///< Loads within the model std::vector loadzip; ///< Loads within the model std::vector gov; ///< Governors within the model + std::vector esdc1a; ///< ESDC1A exciters within the model std::vector exciter; ///< Exciters within the model std::vector sexspti; ///< SEXS-PTI exciters within the model std::vector stabilizer; ///< Stabilizers within the model diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index 00a6278c2..0b882c9df 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -148,6 +148,12 @@ namespace GridKit raw_component.get_to(exciter); sm.exciter.push_back(exciter); } + else if (kind == "Esdc1a") + { + typename SystemModelData::Esdc1aDataT exciter; + raw_component.get_to(exciter); + sm.esdc1a.push_back(exciter); + } else if (kind == "SexsPti") { typename SystemModelData::SexsPtiDataT exciter; diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 932cb13b8..41f8623a1 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -86,6 +86,12 @@ target_link_libraries( GridKit::phasor_dynamics_components_dependency_tracking GridKit::testing) +add_executable(test_phasor_exciter_esdc1a runExciterEsdc1aTests.cpp) +target_link_libraries(test_phasor_exciter_esdc1a GridKit::definitions + GridKit::phasor_dynamics_components + GridKit::phasor_dynamics_components_dependency_tracking + GridKit::testing) + add_executable(test_phasor_exciter_sexspti runExciterSexsPtiTests.cpp) target_link_libraries( test_phasor_exciter_sexspti @@ -136,6 +142,7 @@ add_test(NAME PhasorDynamicsBranchTest COMMAND test_phasor_branch) add_test(NAME PhasorDynamicsGenrouTest COMMAND test_phasor_genrou) add_test(NAME PhasorDynamicsGovernorTgov1Test COMMAND test_phasor_governor_tgov1) add_test(NAME PhasorDynamicsExciterIeeet1Test COMMAND test_phasor_exciter_ieeet1) +add_test(NAME PhasorDynamicsExciterEsdc1aTest COMMAND test_phasor_exciter_esdc1a) add_test(NAME PhasorDynamicsGensalTest COMMAND test_phasor_gensal) add_test(NAME PhasorDynamicsExciterSexsPtiTest COMMAND test_phasor_exciter_sexspti) add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) @@ -157,6 +164,7 @@ install( test_phasor_genrou test_phasor_governor_tgov1 test_phasor_exciter_ieeet1 + test_phasor_exciter_esdc1a test_phasor_gensal test_phasor_exciter_sexspti test_phasor_stabilizer_ieeest diff --git a/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp b/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp new file mode 100644 index 000000000..a2d6a5040 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp @@ -0,0 +1,585 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class ExciterEsdc1aTests + { + public: + using RealT = typename PhasorDynamics::Component::RealT; + + ExciterEsdc1aTests() = default; + ~ExciterEsdc1aTests() = default; + + static constexpr ScalarT kTol = static_cast(1.0e-12); + + TestOutcome constructor() + { + using namespace PhasorDynamics::Exciter; + + TestStatus success = true; + + PhasorDynamics::Bus bus(3.0, 4.0); + + Esdc1a default_exciter(&bus); + success *= (default_exciter.size() == static_cast(Esdc1aInternalVariables::MAXIMUM)); + success *= (default_exciter.getMonitor() == nullptr); + + auto data = makeDefaultData(); + data.monitored_variables.insert(Esdc1aMonitorableVariables::efd); + Esdc1a data_exciter(&bus, data); + success *= (data_exciter.size() == static_cast(Esdc1aInternalVariables::MAXIMUM)); + success *= (data_exciter.getMonitor() != nullptr); + + PhasorDynamics::SignalNode efd_node; + ScalarT efd_value{0.0}; + IdxT efd_index = INVALID_INDEX; + efd_node.set(&efd_value, &efd_index); + + data_exciter.getSignals().template assignSignalNode(&efd_node); + data_exciter.allocate(); + data_exciter.tagDifferentiable(); + + success *= (data_exciter.verify() == 0); + success *= (data_exciter.tag()[idx(Esdc1aInternalVariables::EFDP)] == true); + success *= (data_exciter.tag()[idx(Esdc1aInternalVariables::VC)] == false); + success *= (data_exciter.tag()[idx(Esdc1aInternalVariables::VR)] == true); + success *= (data_exciter.tag()[idx(Esdc1aInternalVariables::VF)] == true); + success *= (data_exciter.tag()[idx(Esdc1aInternalVariables::XLL)] == false); + + return success.report(__func__); + } + + TestOutcome zeroInitialResidual() + { + using namespace PhasorDynamics::Exciter; + + TestStatus success = true; + + Fixture fixture(makeDefaultData()); + success *= fixture.allocateAndInitialize(1.2); + + const auto& f = fixture.exciter.getResidual(); + for (size_t i = 0; i < f.size(); ++i) + { + if (!isEqual(f[i], static_cast(0.0), kTol)) + { + std::cout << "Non-zero ESDC1A residual at index " << i << ": " << f[i] << "\n"; + success = false; + } + } + + success *= fixture.efd_node.linked(); + success *= (fixture.efd_node.getVariableIndex() + == static_cast(idx(Esdc1aInternalVariables::EFD))); + success *= isEqual(fixture.efd_node.read(), static_cast(1.2), kTol); + + return success.report(__func__); + } + + TestOutcome blockDiagramSemantics() + { + TestStatus success = true; + + success *= voltageErrorSummingJunction(); + success *= speedMultiplierSelector(); + success *= leadLagSelector(); + success *= uelRoutingSelector(); + success *= exciterFeedbackLimiter(); + + return success.report(__func__); + } + + TestOutcome parameterValidation() + { + using namespace PhasorDynamics::Exciter; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + PhasorDynamics::SignalNode efd_node; + ScalarT efd_value{0.0}; + IdxT efd_index = INVALID_INDEX; + efd_node.set(&efd_value, &efd_index); + + auto valid = makeDefaultData(); + Esdc1a valid_model(&bus, valid); + valid_model.getSignals().template assignSignalNode(&efd_node); + valid_model.allocate(); + success *= (valid_model.verify() == 0); + + auto invalid_ta = makeDefaultData(); + invalid_ta.parameters[Esdc1aParameters::Ta] = 0.0; + Esdc1a invalid_ta_model(&bus, invalid_ta); + invalid_ta_model.getSignals().template assignSignalNode(&efd_node); + invalid_ta_model.allocate(); + success *= (invalid_ta_model.verify() > 0); + + auto invalid_lead_lag = makeDefaultData(); + invalid_lead_lag.parameters[Esdc1aParameters::Tc] = 0.1; + Esdc1a invalid_lead_lag_model(&bus, invalid_lead_lag); + invalid_lead_lag_model.getSignals().template assignSignalNode(&efd_node); + invalid_lead_lag_model.allocate(); + success *= (invalid_lead_lag_model.verify() > 0); + + auto missing_efd = makeDefaultData(); + Esdc1a missing_efd_model(&bus, missing_efd); + missing_efd_model.allocate(); + success *= (missing_efd_model.verify() > 0); + + auto missing_speed = makeDefaultData(); + missing_speed.parameters[Esdc1aParameters::Spdmlt] = 1.0; + Esdc1a missing_speed_model(&bus, missing_speed); + missing_speed_model.getSignals().template assignSignalNode(&efd_node); + missing_speed_model.allocate(); + success *= (missing_speed_model.verify() > 0); + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + TestOutcome jacobianStructureAndValues() + { + TestStatus success = true; + + const auto tol = static_cast(1.0e-9); + + auto data = makeDefaultData(); + data.parameters[Params::Spdmlt] = 1.0; + data.parameters[Params::UEL] = static_cast(2); + + auto dependency_tracking_jacobian = dependencyTrackingJacobian(data); + auto enzyme_jacobian = enzymeJacobian(data); + + success *= (dependency_tracking_jacobian.size() == enzyme_jacobian.size()); + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) + { + success *= isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], tol); + } + + return success.report(__func__); + } +#endif + + private: + using Internal = PhasorDynamics::Exciter::Esdc1aInternalVariables; + using External = PhasorDynamics::Exciter::Esdc1aExternalVariables; + using Params = PhasorDynamics::Exciter::Esdc1aParameters; + using DataT = PhasorDynamics::Exciter::Esdc1aData; + + static size_t idx(Internal variable) + { + return static_cast(variable); + } + + auto makeDefaultData() -> DataT + { + DataT data; + data.device_class = "exciter"; + data.disambiguation_string = "esdc1a_test"; + + data.parameters[Params::Tr] = 0.0; + data.parameters[Params::Ka] = 40.0; + data.parameters[Params::Ta] = 0.1; + data.parameters[Params::Tb] = 0.0; + data.parameters[Params::Tc] = 0.0; + data.parameters[Params::Vrmax] = 1.0; + data.parameters[Params::Vrmin] = -1.0; + data.parameters[Params::Ke] = 0.1; + data.parameters[Params::Te] = 0.5; + data.parameters[Params::Kf] = 0.05; + data.parameters[Params::Tf1] = 0.7; + data.parameters[Params::Spdmlt] = 0.0; + data.parameters[Params::E1] = 2.8; + data.parameters[Params::Se1] = 0.08; + data.parameters[Params::E2] = 3.7; + data.parameters[Params::Se2] = 0.33; + data.parameters[Params::UEL] = static_cast(0); + data.parameters[Params::exclim] = 1.0; + + return data; + } + + struct Fixture + { + using BusT = PhasorDynamics::Bus; + using SignalT = PhasorDynamics::SignalNode; + using ExciterT = PhasorDynamics::Exciter::Esdc1a; + + DataT data; + BusT bus; + SignalT efd_node; + SignalT omega_node; + SignalT vs_node; + SignalT vuel_node; + + ScalarT efd_value{0.0}; + ScalarT omega_value{0.0}; + ScalarT vs_value{0.0}; + ScalarT vuel_value{-2.0}; + + IdxT efd_index{INVALID_INDEX}; + IdxT omega_index{20}; + IdxT vs_index{21}; + IdxT vuel_index{22}; + + ExciterT exciter; + + explicit Fixture(const DataT& data_in) + : data(data_in), + bus(3.0, 4.0), + exciter(&bus, data) + { + efd_node.set(&efd_value, &efd_index); + omega_node.set(&omega_value, &omega_index); + vs_node.set(&vs_value, &vs_index); + vuel_node.set(&vuel_value, &vuel_index); + + exciter.getSignals().template assignSignalNode(&efd_node); + exciter.getSignals().template attachSignalNode(&omega_node); + exciter.getSignals().template attachSignalNode(&vs_node); + exciter.getSignals().template attachSignalNode(&vuel_node); + } + + bool allocateAndInitialize(ScalarT efd0) + { + bus.allocate(); + bus.initialize(); + exciter.allocate(); + efd_node.init(efd0); + return exciter.verify() == 0 + && exciter.initialize() == 0 + && exciter.evaluateResidual() == 0; + } + }; + + bool voltageErrorSummingJunction() + { + Fixture fixture(makeDefaultData()); + if (!fixture.allocateAndInitialize(1.2)) + { + return false; + } + + fixture.vs_value += 0.1; + fixture.exciter.evaluateResidual(); + bool success = fixture.exciter.getResidual()[idx(Internal::EV)] > static_cast(0.0); + + fixture.vs_value -= 0.1; + fixture.exciter.y()[idx(Internal::VC)] += 0.1; + fixture.exciter.evaluateResidual(); + success = success && fixture.exciter.getResidual()[idx(Internal::EV)] < static_cast(0.0); + + fixture.exciter.y()[idx(Internal::VC)] -= 0.1; + fixture.exciter.y()[idx(Internal::VF)] += 0.1; + fixture.exciter.evaluateResidual(); + success = success && fixture.exciter.getResidual()[idx(Internal::EV)] < static_cast(0.0); + + return success; + } + + bool speedMultiplierSelector() + { + auto disabled_data = makeDefaultData(); + Fixture disabled(disabled_data); + if (!disabled.allocateAndInitialize(1.2)) + { + return false; + } + + disabled.omega_value = 0.05; + disabled.exciter.evaluateResidual(); + bool success = isEqual(disabled.exciter.getResidual()[idx(Internal::EFD)], + static_cast(0.0), + kTol); + + auto enabled_data = makeDefaultData(); + enabled_data.parameters[Params::Spdmlt] = 1.0; + Fixture enabled(enabled_data); + if (!enabled.allocateAndInitialize(1.2)) + { + return false; + } + + enabled.omega_value = 0.05; + enabled.exciter.evaluateResidual(); + success = success && enabled.exciter.getResidual()[idx(Internal::EFD)] > static_cast(0.0); + + return success; + } + + bool leadLagSelector() + { + Fixture bypass(makeDefaultData()); + if (!bypass.allocateAndInitialize(1.2)) + { + return false; + } + + bypass.exciter.y()[idx(Internal::VLL)] += 0.1; + bypass.exciter.evaluateResidual(); + bool success = bypass.exciter.getResidual()[idx(Internal::VLL)] < static_cast(0.0); + + auto active_data = makeDefaultData(); + active_data.parameters[Params::Tb] = 0.5; + active_data.parameters[Params::Tc] = 0.2; + Fixture active(active_data); + if (!active.allocateAndInitialize(1.2)) + { + return false; + } + + active.exciter.y()[idx(Internal::EV)] += 0.1; + active.exciter.evaluateResidual(); + success = success && active.exciter.getResidual()[idx(Internal::VLL)] > static_cast(0.0); + + active.exciter.y()[idx(Internal::EV)] -= 0.1; + active.exciter.y()[idx(Internal::VLL)] += 0.1; + active.exciter.evaluateResidual(); + success = success && active.exciter.getResidual()[idx(Internal::VLL)] < static_cast(0.0); + + return success; + } + + bool uelRoutingSelector() + { + Fixture hv_gate(makeDefaultData()); + if (!hv_gate.allocateAndInitialize(1.2)) + { + return false; + } + + hv_gate.vuel_value = hv_gate.exciter.y()[idx(Internal::VLL)] + 0.1; + hv_gate.exciter.evaluateResidual(); + bool success = hv_gate.exciter.getResidual()[idx(Internal::VHV)] > static_cast(0.0); + success = success && isEqual(hv_gate.exciter.getResidual()[idx(Internal::EV)], static_cast(0.0), kTol); + + auto sum_data = makeDefaultData(); + sum_data.parameters[Params::UEL] = static_cast(2); + Fixture sum_junction(sum_data); + sum_junction.vuel_value = 0.0; + if (!sum_junction.allocateAndInitialize(1.2)) + { + return false; + } + + sum_junction.vuel_value = 0.1; + sum_junction.exciter.evaluateResidual(); + success = success && sum_junction.exciter.getResidual()[idx(Internal::EV)] > static_cast(0.0); + success = success && isEqual(sum_junction.exciter.getResidual()[idx(Internal::VHV)], static_cast(0.0), kTol); + + return success; + } + + bool exciterFeedbackLimiter() + { + auto limited_data = makeDefaultData(); + limited_data.parameters[Params::Ke] = -0.2; + limited_data.parameters[Params::Se1] = 0.0; + limited_data.parameters[Params::Se2] = 0.0; + limited_data.parameters[Params::exclim] = 1.0; + Fixture limited(limited_data); + if (!limited.allocateAndInitialize(1.2)) + { + return false; + } + + limited.exciter.y()[idx(Internal::EFDP)] = 1.0; + limited.exciter.y()[idx(Internal::SE)] = 0.0; + limited.exciter.y()[idx(Internal::VFE)] = 0.0; + limited.exciter.evaluateResidual(); + bool success = std::abs(limited.exciter.getResidual()[idx(Internal::VFE)]) < kTol; + + auto unlimited_data = limited_data; + unlimited_data.parameters[Params::exclim] = 0.0; + Fixture unlimited(unlimited_data); + if (!unlimited.allocateAndInitialize(1.2)) + { + return false; + } + + unlimited.exciter.y()[idx(Internal::EFDP)] = 1.0; + unlimited.exciter.y()[idx(Internal::SE)] = 0.0; + unlimited.exciter.y()[idx(Internal::VFE)] = 0.0; + unlimited.exciter.evaluateResidual(); + success = success && unlimited.exciter.getResidual()[idx(Internal::VFE)] < static_cast(0.0); + + return success; + } + +#ifdef GRIDKIT_ENABLE_ENZYME + std::vector + dependencyTrackingJacobian(const DataT& data) + { + using Variable = DependencyTracking::Variable; + + PhasorDynamics::Bus bus(Variable{3.0}, Variable{4.0}); + PhasorDynamics::SignalNode efd_node; + PhasorDynamics::SignalNode omega_node; + PhasorDynamics::SignalNode vs_node; + PhasorDynamics::SignalNode vuel_node; + + Variable efd_value{0.0}; + Variable omega_value{0.0}; + Variable vs_value{0.0}; + Variable vuel_value{0.0}; + + IdxT efd_index = INVALID_INDEX; + IdxT omega_index = 13; + IdxT vs_index = 14; + IdxT vuel_index = 15; + + efd_node.set(&efd_value, &efd_index); + omega_node.set(&omega_value, &omega_index); + vs_node.set(&vs_value, &vs_index); + vuel_node.set(&vuel_value, &vuel_index); + + PhasorDynamics::Exciter::Esdc1a exciter(&bus, data); + exciter.getSignals().template assignSignalNode(&efd_node); + exciter.getSignals().template attachSignalNode(&omega_node); + exciter.getSignals().template attachSignalNode(&vs_node); + exciter.getSignals().template attachSignalNode(&vuel_node); + + bus.allocate(); + exciter.allocate(); + bus.initialize(); + efd_node.init(Variable{1.2}); + exciter.initialize(); + + for (size_t i = 0; i < exciter.size(); ++i) + { + exciter.y()[i].setVariableNumber(i); + } + for (size_t i = 0; i < bus.size(); ++i) + { + bus.y()[i].setVariableNumber(i + exciter.size()); + } + omega_value.setVariableNumber(13); + vs_value.setVariableNumber(14); + vuel_value.setVariableNumber(15); + + bus.evaluateResidual(); + exciter.evaluateResidual(); + auto residual_y = exciter.getResidual(); + + omega_value = 0.0; + vs_value = 0.0; + vuel_value = 0.0; + bus.initialize(); + efd_node.init(Variable{1.2}); + exciter.initialize(); + + for (size_t i = 0; i < exciter.size(); ++i) + { + exciter.yp()[i].setVariableNumber(i); + } + + bus.evaluateResidual(); + exciter.evaluateResidual(); + auto residual_yp = exciter.getResidual(); + + std::vector dependencies(residual_y.size()); + for (IdxT i = 0; i < residual_y.size(); ++i) + { + auto dependency_y = residual_y[static_cast(i)].getDependencies(); + auto dependency_yp = residual_yp[static_cast(i)].getDependencies(); + + for (const auto& pair_y : dependency_y) + { + auto index_y = pair_y.first; + auto value_y = pair_y.second; + auto it_yp = dependency_yp.find(index_y); + if (it_yp != dependency_yp.end()) + { + dependencies[static_cast(i)].insert(std::make_pair(index_y, value_y + it_yp->second)); + } + else + { + dependencies[static_cast(i)].insert(pair_y); + } + } + + for (const auto& pair_yp : dependency_yp) + { + if (dependency_y.find(pair_yp.first) == dependency_y.end()) + { + dependencies[static_cast(i)].insert(pair_yp); + } + } + } + + return dependencies; + } + + std::vector + enzymeJacobian(const DataT& data) + { + PhasorDynamics::Bus bus(3.0, 4.0); + PhasorDynamics::SignalNode efd_node; + PhasorDynamics::SignalNode omega_node; + PhasorDynamics::SignalNode vs_node; + PhasorDynamics::SignalNode vuel_node; + + ScalarT efd_value{0.0}; + ScalarT omega_value{0.0}; + ScalarT vs_value{0.0}; + ScalarT vuel_value{0.0}; + + IdxT efd_index = INVALID_INDEX; + IdxT omega_index = 13; + IdxT vs_index = 14; + IdxT vuel_index = 15; + + efd_node.set(&efd_value, &efd_index); + omega_node.set(&omega_value, &omega_index); + vs_node.set(&vs_value, &vs_index); + vuel_node.set(&vuel_value, &vuel_index); + + PhasorDynamics::Exciter::Esdc1a exciter(&bus, data); + exciter.getSignals().template assignSignalNode(&efd_node); + exciter.getSignals().template attachSignalNode(&omega_node); + exciter.getSignals().template attachSignalNode(&vs_node); + exciter.getSignals().template attachSignalNode(&vuel_node); + + bus.allocate(); + exciter.allocate(); + bus.initialize(); + efd_node.init(1.2); + exciter.initialize(); + exciter.updateTime(0.0, 1.0); + + for (size_t i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, static_cast(i + exciter.size())); + bus.setResidualIndex(i, static_cast(i + exciter.size())); + } + + bus.evaluateResidual(); + exciter.evaluateResidual(); + exciter.evaluateJacobian(); + + auto& model_jacobian = exciter.getJacobian(); + model_jacobian.deduplicate(); + + return MapFromCOO(model_jacobian); + } +#endif + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runExciterEsdc1aTests.cpp b/tests/UnitTests/PhasorDynamics/runExciterEsdc1aTests.cpp new file mode 100644 index 000000000..04c40b955 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runExciterEsdc1aTests.cpp @@ -0,0 +1,18 @@ +#include "ExciterEsdc1aTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::ExciterEsdc1aTests test; + + result += test.constructor(); + result += test.zeroInitialResidual(); + result += test.blockDiagramSemantics(); + result += test.parameterValidation(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobianStructureAndValues(); +#endif + + return result.summary(); +} From 7867706f8173fde8ca10b36ecf1a54238a004fd3 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Tue, 2 Jun 2026 20:16:40 +0000 Subject: [PATCH 3/3] Apply pre-commmit fixes --- .../Exciter/ESDC1A/CMakeLists.txt | 65 +++++++++---------- tests/UnitTests/PhasorDynamics/CMakeLists.txt | 10 +-- 2 files changed, 37 insertions(+), 38 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt index 6886a26ec..53807c5df 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt @@ -3,47 +3,44 @@ # - Luke Lowery # ]] -set(_install_headers - Esdc1a.hpp - Esdc1aData.hpp) +set(_install_headers Esdc1a.hpp Esdc1aData.hpp) if(GRIDKIT_ENABLE_ENZYME) - gridkit_add_library(phasor_dynamics_exciter_esdc1a - SOURCES - Esdc1aEnzyme.cpp - HEADERS - ${_install_headers} - INCLUDE_DIRECTORIES - PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + gridkit_add_library( + phasor_dynamics_exciter_esdc1a + SOURCES Esdc1aEnzyme.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include LINK_LIBRARIES - PUBLIC GridKit::phasor_dynamics_core - PUBLIC GridKit::phasor_dynamics_signal - PRIVATE ClangEnzymeFlags + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal + PRIVATE + ClangEnzymeFlags COMPILE_OPTIONS - PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno) + PRIVATE + -mllvm + -enzyme-auto-sparsity=1 + -fno-math-errno) else() - gridkit_add_library(phasor_dynamics_exciter_esdc1a - SOURCES - Esdc1a.cpp - HEADERS - ${_install_headers} - INCLUDE_DIRECTORIES - PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include - LINK_LIBRARIES - GridKit::phasor_dynamics_core - GridKit::phasor_dynamics_signal) + gridkit_add_library( + phasor_dynamics_exciter_esdc1a + SOURCES Esdc1a.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES GridKit::phasor_dynamics_core GridKit::phasor_dynamics_signal) endif() -gridkit_add_library(phasor_dynamics_exciter_esdc1a_dependency_tracking - SOURCES - Esdc1aDependencyTracking.cpp - INCLUDE_DIRECTORIES - PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include - LINK_LIBRARIES - GridKit::phasor_dynamics_core - GridKit::phasor_dynamics_signal_dependency_tracking) +gridkit_add_library( + phasor_dynamics_exciter_esdc1a_dependency_tracking + SOURCES Esdc1aDependencyTracking.cpp + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES GridKit::phasor_dynamics_core GridKit::phasor_dynamics_signal_dependency_tracking) -target_link_libraries(phasor_dynamics_components +target_link_libraries( + phasor_dynamics_components INTERFACE GridKit::phasor_dynamics_exciter_esdc1a) -target_link_libraries(phasor_dynamics_components_dependency_tracking +target_link_libraries( + phasor_dynamics_components_dependency_tracking INTERFACE GridKit::phasor_dynamics_exciter_esdc1a_dependency_tracking) diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 41f8623a1..451f18966 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -87,10 +87,12 @@ target_link_libraries( GridKit::testing) add_executable(test_phasor_exciter_esdc1a runExciterEsdc1aTests.cpp) -target_link_libraries(test_phasor_exciter_esdc1a GridKit::definitions - GridKit::phasor_dynamics_components - GridKit::phasor_dynamics_components_dependency_tracking - GridKit::testing) +target_link_libraries( + test_phasor_exciter_esdc1a + GridKit::definitions + GridKit::phasor_dynamics_components + GridKit::phasor_dynamics_components_dependency_tracking + GridKit::testing) add_executable(test_phasor_exciter_sexspti runExciterSexsPtiTests.cpp) target_link_libraries(