From d3b5364709ea2a68971e889d89adf91ec623b052 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Fri, 29 May 2026 23:36:18 -0500 Subject: [PATCH] Add REECB phasor dynamics model --- CHANGELOG.md | 1 + GridKit/CommonMath.hpp | 24 +- .../Model/PhasorDynamics/ComponentLibrary.hpp | 1 + .../PhasorDynamics/Converter/CMakeLists.txt | 1 + .../Model/PhasorDynamics/Converter/README.md | 1 + .../Converter/REECB/CMakeLists.txt | 54 ++ .../PhasorDynamics/Converter/REECB/README.md | 376 +++++++++ .../PhasorDynamics/Converter/REECB/Reecb.cpp | 27 + .../PhasorDynamics/Converter/REECB/Reecb.hpp | 188 +++++ .../Converter/REECB/ReecbData.hpp | 100 +++ .../REECB/ReecbDependencyTracking.cpp | 27 + .../Converter/REECB/ReecbEnzyme.cpp | 92 +++ .../Converter/REECB/ReecbImpl.hpp | 721 ++++++++++++++++++ GridKit/Model/PhasorDynamics/INPUT_FORMAT.md | 4 +- GridKit/Model/PhasorDynamics/SystemModel.hpp | 61 +- .../Model/PhasorDynamics/SystemModelData.hpp | 3 + .../SystemModelDataJSONParser.hpp | 6 + docs/Figures/PhasorDynamics_REECB_Diagram.png | Bin 0 -> 77117 bytes .../Medium/NewEngland/CMakeLists.txt | 1 - tests/UnitTests/PhasorDynamics/CMakeLists.txt | 10 + .../PhasorDynamics/ConverterReecbTests.hpp | 535 +++++++++++++ .../PhasorDynamics/runConverterReecbTests.cpp | 22 + 22 files changed, 2242 insertions(+), 13 deletions(-) create mode 100644 GridKit/Model/PhasorDynamics/Converter/REECB/CMakeLists.txt create mode 100644 GridKit/Model/PhasorDynamics/Converter/REECB/README.md create mode 100644 GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.cpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.hpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REECB/ReecbData.hpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REECB/ReecbDependencyTracking.cpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REECB/ReecbEnzyme.cpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REECB/ReecbImpl.hpp create mode 100644 docs/Figures/PhasorDynamics_REECB_Diagram.png create mode 100644 tests/UnitTests/PhasorDynamics/ConverterReecbTests.hpp create mode 100644 tests/UnitTests/PhasorDynamics/runConverterReecbTests.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index 71751c60a..63963c6d4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -62,6 +62,7 @@ - Added cmake-format hooks, including in pre-commit. - Added off-nominal tap ratio and phase shift support to the PhasorDynamics `Branch` model. - Added `REGCA` converter model implementation for PhasorDynamics. +- Added `REECB` converter model implementation for PhasorDynamics. ## v0.1 diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index daaa11e94..130008eb0 100644 --- a/GridKit/CommonMath.hpp +++ b/GridKit/CommonMath.hpp @@ -15,7 +15,7 @@ namespace GridKit * * @note The sigmoid constant (mu) value is chosen to balance accuracy * and finite derivatives. Large values more closely approximate a step - * function, but lead to inf or NaN derivatives. + * function, but can make the transition numerically stiff. * * @tparam ScalarT - scalar data type * @@ -27,7 +27,7 @@ namespace GridKit { using RealT = typename GridKit::ScalarTraits::RealT; static constexpr RealT MU = 240.0; - return ONE / (ONE + std::exp(-MU * x)); + return HALF * (ONE + std::tanh(HALF * MU * x)); } /** @@ -159,7 +159,8 @@ namespace GridKit * passes the input through unchanged outside the band. * * @tparam ScalarT - scalar data type - * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * @tparam LowerT - data type of the lower limit + * @tparam UpperT - data type of the upper limit * * @param[in] x - Input signal * @param[in] lower - Lower breakpoint @@ -355,13 +356,15 @@ namespace GridKit * @return Scalar value in [0, 1]: 1 when dynamics should pass through, * 0 when integration should be blocked. */ - template + template __attribute__((always_inline)) inline ScalarT indicator( const ScalarT x, const ScalarT f, - const RealT limit_min, - const RealT limit_max) + const LowerT limit_min, + const UpperT limit_max) { + using RealT = typename GridKit::ScalarTraits::RealT; + assert(limit_min <= limit_max); ScalarT above_min = above(x, limit_min); @@ -381,7 +384,8 @@ namespace GridKit * and blocks motion that would push further into saturation. * * @tparam ScalarT - Scalar data type - * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * @tparam LowerT - data type of the lower limit + * @tparam UpperT - data type of the upper limit * * @param[in] x - Limited state or limited output signal * @param[in] f - Pre-limit derivative @@ -389,12 +393,12 @@ namespace GridKit * @param[in] limit_max - Maximum limit * @return Smooth anti-windup limited derivative */ - template + template __attribute__((always_inline)) inline ScalarT antiwindup( const ScalarT x, const ScalarT f, - const RealT limit_min, - const RealT limit_max) + const LowerT limit_min, + const UpperT limit_max) { return indicator(x, f, limit_min, limit_max) * f; } diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index e30067d20..fde8eaa8a 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/Converter/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt index cafb2cc36..5a7dffab8 100644 --- a/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt @@ -4,3 +4,4 @@ # ]] add_subdirectory(REGCA) +add_subdirectory(REECB) diff --git a/GridKit/Model/PhasorDynamics/Converter/README.md b/GridKit/Model/PhasorDynamics/Converter/README.md index ad38ba19a..30b760ecb 100644 --- a/GridKit/Model/PhasorDynamics/Converter/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/README.md @@ -12,3 +12,4 @@ The GridKit converter documentation includes: - Renewable Energy Generator/Converter Model REGCA (See [REGCA](REGCA/README.md)) - Renewable Energy Generator/Converter Model REGCB (See [REGCB](REGCB/README.md)) - Renewable Energy Electrical Control Model REECA (See [REECA](REECA/README.md)) +- Renewable Energy Electrical Control Model REECB (See [REECB](REECB/README.md)) diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/REECB/CMakeLists.txt new file mode 100644 index 000000000..03492a503 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Reecb.hpp ReecbData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_converter_reecb + SOURCES ReecbEnzyme.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_converter_reecb + SOURCES Reecb.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) +endif() + +gridkit_add_library( + phasor_dynamics_converter_reecb_dependency_tracking + SOURCES ReecbDependencyTracking.cpp + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal_dependency_tracking) + +target_link_libraries( + phasor_dynamics_components + INTERFACE GridKit::phasor_dynamics_converter_reecb) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_converter_reecb_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/README.md b/GridKit/Model/PhasorDynamics/Converter/REECB/README.md new file mode 100644 index 000000000..fa8ce3e33 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/README.md @@ -0,0 +1,376 @@ +# **Renewable Energy Electrical Control Model (REECB)** + +REECB is a WECC renewable energy electrical control model for inverter-coupled +resources. In GridKit it is represented as a signal-control model that computes +active- and reactive-current commands. + +Notes: +- Internal electrical quantities and current commands are on model base unless otherwise stated. +- `pe`/`qgen` electrical signal ports are on system base and are converted to + REECB model base through `mva`. +- The `pe` and `qgen` signal ports must be connected together. If both are + omitted, initialization recovers the constant electrical feedback from the + model-base `ipcmd` and `iqcmd` signals. +- `qext`, `pref`, `iqcmd`, and `ipcmd` are model-base command signals. +- REECB uses the REECA voltage-dip freeze structure, but does not include the REECA VDL current-limit curves, active-power speed multiplier, or post-dip hold timers. +- Source governor-response settings may modify active-power ramp-rate or order limits before the model equations are evaluated; the equations below use those effective limits. + +## Block Diagram + +Standard REECB block diagram. + +
+ + + Figure 1: REECB block diagram. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +
+ +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value | Note +------------------------------------|----------|----------|---------------------------------------------------------|---------------|------ +$S^{\mathrm{base}}$ | [MVA] | `mva` | REECB model power base | 0.0 | Block name: `MVABase` +$s_{\mathrm{pf}}$ | [binary] | `PfFlag` | Power-factor control flag | 0 | Block name: `PfFlag`; 1 = power-factor control, 0 = Q control +$s_V$ | [binary] | `VFlag` | Voltage-control mode flag | 0 | Block name: `VFlag`; 1 = Q control, 0 = voltage control +$s_Q$ | [binary] | `QFlag` | Reactive-power control flag | 0 | Block name: `QFlag`; 1 = voltage/Q control, 0 = constant pf or Q control +$s_{PQ}$ | [binary] | `Pqflag` | P/Q priority flag for converter current limit | 0 | Block name: `Pqflag`; 0 = Q priority, 1 = P priority +$T_{\mathrm{rv}}$ | [sec] | `Trv` | Voltage-measurement filter time constant | 0.02 | Block name: `Trv`; if zero, $V_{\mathrm{meas}}$ is algebraic +$T_{\mathrm{p}}$ | [sec] | `Tp` | Electrical-power measurement filter time constant | 0.0 | Block name: `Tp`; if zero, $P_{\mathrm{meas}}$ is algebraic +$V_{\mathrm{ref0}}$ | [p.u.] | `Vref0` | Outer-loop voltage reference | 0.0 | Block name: `Vref0`; initialized to terminal voltage if omitted +$V_{\mathrm{dip}}$ | [p.u.] | `Vdip` | Low-voltage threshold for reactive-current injection logic | 0.85 | Block name: `Vdip` +$V_{\mathrm{up}}$ | [p.u.] | `Vup` | High-voltage threshold for reactive-current injection logic | 1.15 | Block name: `Vup` +$D_{\mathrm{bd1}}$ | [p.u.] | `dbd1` | Overvoltage deadband for voltage-error response | 0.0 | Block name: `dbd1` +$D_{\mathrm{bd2}}$ | [p.u.] | `dbd2` | Undervoltage deadband for voltage-error response | 0.0 | Block name: `dbd2` +$K_{\mathrm{qv}}$ | [p.u.] | `kqv` | Reactive-current injection gain during voltage dip/overvoltage logic | 5.0 | Block name: `kqv` +$I_{\mathrm{qinj}}^{\min}$ | [p.u.] | `Iql1` | Minimum reactive-current injection limit | -1.1 | Block name: `Iql1` +$I_{\mathrm{qinj}}^{\max}$ | [p.u.] | `Iqh1` | Maximum reactive-current injection limit | 1.1 | Block name: `Iqh1` +$Q^{\max}$ | [p.u.] | `Qmax` | Maximum reactive-power control limit | 0.436 | Block name: `Qmax` +$Q^{\min}$ | [p.u.] | `Qmin` | Minimum reactive-power control limit | -0.436 | Block name: `Qmin` +$K_{\mathrm{qp}}$ | [p.u.] | `Kqp` | Reactive-power control proportional gain | 0.0 | Block name: `Kqp` +$K_{\mathrm{qi}}$ | [p.u./s] | `Kqi` | Reactive-power control integral gain | 0.1 | Block name: `Kqi` +$V^{\max}$ | [p.u.] | `Vmax` | Maximum voltage-control limit | 1.1 | Block name: `Vmax` +$V^{\min}$ | [p.u.] | `Vmin` | Minimum voltage-control limit | 0.9 | Block name: `Vmin` +$K_{\mathrm{vp}}$ | [p.u.] | `Kvp` | Voltage-control proportional gain | 18.0 | Block name: `Kvp` +$K_{\mathrm{vi}}$ | [p.u./s] | `Kvi` | Voltage-control integral gain | 5.0 | Block name: `Kvi` +$T_{\mathrm{iq}}$ | [sec] | `Tiq` | Reactive-current command lag time constant | 0.02 | Block name: `Tiq` +$T_{\mathrm{pord}}$ | [sec] | `Tpord` | Active-power order filter time constant | 0.0 | Block name: `Tpord` +$R_P^{\max}$ | [p.u./s] | `dPmax` | Positive active-power order ramp-rate limit | 99.0 | Block name: `dPmax` +$R_P^{\min}$ | [p.u./s] | `dPmin` | Negative active-power order ramp-rate limit | -99.0 | Block name: `dPmin` +$P^{\max}$ | [p.u.] | `Pmax` | Maximum active-power order limit | 1.0 | Block name: `Pmax` +$P^{\min}$ | [p.u.] | `Pmin` | Minimum active-power order limit | 0.0 | Block name: `Pmin` +$I^{\max}$ | [p.u.] | `Imax` | Maximum total converter current | 1.3 | Block name: `Imax` + +### Parameter Validation + +Invalid REECB parameter sets are rejected by the following checks. If source data preprocessing adjusts active-power ramp-rate or order limits, apply these checks to the effective values used by the equations. + +The required checks are: + +```math +\begin{aligned} + &S^{\mathrm{base}} > 0 \\ + &s_{\mathrm{pf}}, s_V, s_Q, s_{PQ} \in \{0,1\} \\ + &T_{\mathrm{rv}}, T_{\mathrm{p}} \ge 0 \\ + &V_{\mathrm{dip}} < V_{\mathrm{up}} \\ + &D_{\mathrm{bd1}} \le 0 \le D_{\mathrm{bd2}} \\ + &I_{\mathrm{qinj}}^{\min} \le I_{\mathrm{qinj}}^{\max} \\ + &Q^{\min} \le Q^{\max} \\ + &V^{\min} \le V^{\max} \\ + &T_{\mathrm{iq}}, T_{\mathrm{pord}} > 0 \\ + &R_P^{\min} < 0 < R_P^{\max} \\ + &P^{\min} \le P^{\max} \\ + &I^{\max} \ge 0 +\end{aligned} +``` + +### Model Derived Parameters + +The off-mode flag complements are: + +```math +\begin{aligned} + s_{\mathrm{pf}}^{\mathrm{off}} &= 1 - s_{\mathrm{pf}} \\ + s_V^{\mathrm{off}} &= 1 - s_V \\ + s_Q^{\mathrm{off}} &= 1 - s_Q \\ + s_{PQ}^{\mathrm{off}} &= 1 - s_{PQ} +\end{aligned} +``` + +System-base electrical signals are converted to REECB model base by: + +```math +\begin{aligned} + x_{\mathrm{model}} &= x_{\mathrm{system}} + \dfrac{S^{\mathrm{sys}}}{S^{\mathrm{base}}} +\end{aligned} +``` + +The same conversion is applied to connected `pe`/`qgen` signals. + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------|--------|-------------------------------------|------ +$V_{\mathrm{meas}}$ | [p.u.] | Filtered terminal voltage | State 1 in Fig. 1; source label: `Vmeas`; algebraic when $T_{\mathrm{rv}} = 0$ +$P_{\mathrm{meas}}$ | [p.u.] | Filtered electrical power | State 2 in Fig. 1; source label: `Pmeas`; algebraic when $T_{\mathrm{p}} = 0$ +$x_{\mathrm{PIQ}}$ | [p.u.] | Reactive-power PI controller state | State 3 in Fig. 1; source label: `PIQ` +$x_{\mathrm{PIV}}$ | [p.u.] | Voltage PI controller state | State 4 in Fig. 1; source label: `PIV` +$Q_V$ | [p.u.] | Reactive-current command lag state | State 5 in Fig. 1; source label: `Q_V` +$P_{\mathrm{ord}}$ | [p.u.] | Filtered active-power order | State 6 in Fig. 1; source label: `Pord` + +#### Algebraic + +Symbol | Units | Description | Note +--------------------------------|--------|-------------------------------------|------ +$V_T$ | [p.u.] | Terminal voltage magnitude | +$V_{\mathrm{meas}}^{\mathrm{safe}}$ | [p.u.] | Safe filtered terminal voltage for divider blocks | Lower bounded by 0.01 +$s_{\mathrm{dip}}$ | [binary] | Voltage-dip/overvoltage freeze indicator | 1 when outside voltage thresholds +$V_{\mathrm{err}}$ | [p.u.] | Deadbanded voltage error | Defined by CommonMath `deadband2` +$I_{\mathrm{qv}}$ | [p.u.] | Reactive-current injection candidate | Converter base +$Q_{\mathrm{ref}}$ | [p.u.] | Selected reactive-power reference | From power-factor or external reactive-power command +$e_Q$ | [p.u.] | Reactive-power control error | Limited $Q_{\mathrm{ref}}$ minus $Q_{\mathrm{gen}}$ +$V_{\mathrm{PIQ}}$ | [p.u.] | Reactive-power control PI output | Limited by $V^{\min}$ and $V^{\max}$ +$e_{\mathrm{PIV}}$ | [p.u.] | Voltage-control PI error | Selected voltage-control signal minus $V_{\mathrm{meas}}$ +$f_{\mathrm{pord}}$ | [p.u./s] | Active-power order derivative before ramp-rate limiting | Feeds $r_{\mathrm{pord}}$ +$r_{\mathrm{pord}}$ | [p.u./s] | Ramp-rate-limited active-power order derivative | Feeds $P_{\mathrm{ord}}$ anti-windup +$I_{\mathrm{q}}^{\mathrm{circ}}$ | [p.u.] | Reactive-current limit from converter current circle | Converter base; nonnegative algebraic branch +$I_{\mathrm{p}}^{\mathrm{circ}}$ | [p.u.] | Active-current limit from converter current circle | Converter base; nonnegative algebraic branch +$I_{\mathrm{q}}^{\max}$ | [p.u.] | Final reactive-current upper limit | Converter base; updated by current-limit logic +$I_{\mathrm{p}}^{\max}$ | [p.u.] | Final active-current upper limit | Converter base; updated by current-limit logic +$I_{\mathrm{qbase}}$ | [p.u.] | Base reactive-current command | Converter base; before $s_Q$ selection and reactive-current injection +$I_{\mathrm{q}}^{\mathrm{raw}}$ | [p.u.] | Raw reactive-current command before final limit | Converter base +$I_{\mathrm{q}}^{\mathrm{cmd}}$ | [p.u.] | Reactive-current command output | Converter base +$I_{\mathrm{p}}^{\mathrm{cmd}}$ | [p.u.] | Active-current command output | Converter base + +### External Variables + +#### Differential +None. + +#### Algebraic + +Symbol | Units | Description | Note +-------------------------------------|--------|-----------------------------------|------ +$V_{\mathrm{r}}$ | [p.u.] | Terminal voltage, real component | Owned by bus object +$V_{\mathrm{i}}$ | [p.u.] | Terminal voltage, imaginary component | Owned by bus object +$P_e$ | [p.u.] | Electrical active power | Source label: `Pe`; system-base signal converted to model base +$Q_{\mathrm{gen}}$ | [p.u.] | Reactive-power signal | Source label: `Qgen`; system-base signal converted to model base +$Q_{\mathrm{ext}}$ | [p.u.] | External reactive-power command | Model-base command; optional, defaults to initialized constant; source label: `Qext` +$\phi_{\mathrm{pf}}^{\mathrm{ref}}$ | [rad] | Power-factor angle reference | Source label: `pfaref`; used through tangent block +$P_{\mathrm{ref}}$ | [p.u.] | External active-power reference | Model-base command; optional, defaults to initialized constant + +## Model Equations + +### Differential Equations + +The state-equation residuals use compact limiter notation where applicable. The measurement filters are written in descriptor form; when $T_{\mathrm{rv}} = 0$ or $T_{\mathrm{p}} = 0$, the corresponding residual becomes algebraic. The reactive-current lag requires $T_{\mathrm{iq}} > 0$. + +```math +\begin{aligned} + 0 &= -T_{\mathrm{rv}}\dot V_{\mathrm{meas}} - V_{\mathrm{meas}} + V_T \\ + 0 &= -T_{\mathrm{p}}\dot P_{\mathrm{meas}} - P_{\mathrm{meas}} + P_e \\ + 0 &= + -\dot x_{\mathrm{PIQ}} + + (1 - s_{\mathrm{dip}}) + \text{antiwindup}\!\left( + V_{\mathrm{PIQ}}, + K_{\mathrm{qi}} e_Q, + V^{\min}, + V^{\max} + \right) \\ + 0 &= + -\dot x_{\mathrm{PIV}} + + (1 - s_{\mathrm{dip}}) + \text{antiwindup}\!\left( + I_{\mathrm{qbase}}, + K_{\mathrm{vi}} e_{\mathrm{PIV}}, + -I_{\mathrm{q}}^{\max}, + I_{\mathrm{q}}^{\max} + \right) \\ + 0 &= + -T_{\mathrm{iq}}\dot Q_V + - (1 - s_{\mathrm{dip}})Q_V + + (1 - s_{\mathrm{dip}})Q_{\mathrm{ref}}/V_{\mathrm{meas}}^{\mathrm{safe}} \\ + 0 &= + -\dot P_{\mathrm{ord}} + + (1 - s_{\mathrm{dip}}) + \text{antiwindup}\!\left( + P_{\mathrm{ord}}, + r_{\mathrm{pord}}, + P^{\min}, + P^{\max} + \right) +\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 &= -V_T^2 + V_\mathrm r^2 + V_\mathrm i^2 \\ + 0 &= -V_\mathrm{meas}^\mathrm{safe} + \max(V_\mathrm{meas}, 0.01) \\ + 0 &= -s_\mathrm{dip} + \text{outside}(V_T, V_\mathrm{dip}, V_\mathrm{up}) \\ + 0 &= -V_\mathrm{err} + \text{deadband2}(V_\mathrm{ref0} - V_\mathrm{meas}, D_\mathrm{bd1}, D_\mathrm{bd2}) \\ + 0 &= -I_\mathrm{qv} + \text{clamp}(K_\mathrm{qv} V_\mathrm{err}, I_\mathrm{qinj}^{\min}, I_\mathrm{qinj}^{\max}) \\ + 0 &= -Q_\mathrm{ref} + + s_\mathrm{pf} P_\mathrm{meas}\tan(\phi_\mathrm{pf}^\mathrm{ref}) + + s_\mathrm{pf}^\mathrm{off} Q_\mathrm{ext} \\ + 0 &= -e_Q + \text{clamp}(Q_\mathrm{ref}, Q^{\min}, Q^{\max}) - Q_\mathrm{gen} \\ + 0 &= -V_\mathrm{PIQ} + \text{clamp}(K_\mathrm{qp} e_Q + x_\mathrm{PIQ}, V^{\min}, V^{\max}) \\ + 0 &= -e_\mathrm{PIV} + s_V V_\mathrm{PIQ} + s_V^\mathrm{off}Q_\mathrm{ref} - V_\mathrm{meas} \\ + 0 &= -T_\mathrm{pord} f_\mathrm{pord} + P_\mathrm{ref} - P_\mathrm{ord} \\ + 0 &= -r_\mathrm{pord} + \text{clamp}(f_\mathrm{pord}, R_P^{\min}, R_P^{\max}) +\end{aligned} +``` + +```math +\begin{aligned} + 0 &= -{I_\mathrm{q}^\mathrm{circ}}^2 + (I^{\max})^2 - s_{PQ}(I_\mathrm{p}^\mathrm{cmd})^2 \\ + 0 &= -{I_\mathrm{p}^\mathrm{circ}}^2 + (I^{\max})^2 - s_{PQ}^\mathrm{off}(I_\mathrm{q}^\mathrm{cmd})^2 \\ + 0 &= -I_\mathrm{q}^{\max} + s_{PQ}^{\mathrm{off}} I^{\max} + s_{PQ} I_\mathrm{q}^{\mathrm{circ}} \\ + 0 &= -I_\mathrm{p}^{\max} + s_{PQ} I^{\max} + s_{PQ}^{\mathrm{off}} I_\mathrm{p}^{\mathrm{circ}} \\ + 0 &= -I_\mathrm{qbase} + \text{clamp}(K_\mathrm{vp} e_\mathrm{PIV} + x_\mathrm{PIV}, -I_\mathrm{q}^{\max}, I_\mathrm{q}^{\max}) \\ + 0 &= -I_\mathrm{q}^\mathrm{raw} + s_Q I_\mathrm{qbase} + s_Q^\mathrm{off} Q_V + s_\mathrm{dip} I_\mathrm{qv} \\ + 0 &= -I_\mathrm{q}^\mathrm{cmd} + \text{clamp}(I_\mathrm{q}^\mathrm{raw}, -I_\mathrm{q}^{\max}, I_\mathrm{q}^{\max}) \\ + 0 &= -I_\mathrm{p}^\mathrm{cmd} + \text{clamp}(P_\mathrm{ord}/V_\mathrm{meas}^\mathrm{safe}, 0, I_\mathrm{p}^{\max}) +\end{aligned} +``` + +The $V_T$, $I_{\mathrm{q}}^{\mathrm{circ}}$, and $I_{\mathrm{p}}^{\mathrm{circ}}$ variables use nonnegative branches of squared algebraic residuals. This matches the REECB current-limit pseudo-code: for $s_{PQ}=0$, $I_{\mathrm{q}}^{\max}=I^{\max}$ and $I_{\mathrm{p}}^{\max}$ is reduced by $I_{\mathrm{q}}^{\mathrm{cmd}}$; for $s_{PQ}=1$, $I_{\mathrm{p}}^{\max}=I^{\max}$ and $I_{\mathrm{q}}^{\max}$ is reduced by $I_{\mathrm{p}}^{\mathrm{cmd}}$. A consistent solution satisfies the nonnegative branch and nonnegative radicands. + +CommonMath defines the helper targets and smooth approximations for [max, clamp, deadband2, and outside](../../../../CommonMath.md#derived-functions). + +## Initialization + +Initialization is performed by evaluating the steady-state residuals in dependency order. Let subscript $0$ denote initial values and set all internal derivatives to zero. REECB recovers the initial electrical operating point from the `ipcmd` and `iqcmd` signal values initialized upstream: + +```math +\begin{aligned} + V_{T,0} &= \sqrt{V_{\mathrm{r},0}^2 + V_{\mathrm{i},0}^2} \\ + V_{\mathrm{meas},0} &= V_{T,0} \\ + V_{\mathrm{meas},0}^{\mathrm{safe}} &= \text{max}(V_{\mathrm{meas},0}, 0.01) \\ + P_{e,0} &= I_{\mathrm{p},0}^{\mathrm{cmd}}V_{\mathrm{meas},0}^{\mathrm{safe}} \\ + Q_{\mathrm{gen},0} &= I_{\mathrm{q},0}^{\mathrm{cmd}}V_{\mathrm{meas},0}^{\mathrm{safe}} \\ + Q_{\mathrm{ext},0} &= Q_{\mathrm{gen},0} \\ + \phi_{\mathrm{pf},0}^{\mathrm{ref}} &= \begin{cases} + \tan^{-1}(Q_{\mathrm{ext},0}/P_{e,0}), & P_{e,0} \neq 0 \\ + 0, & P_{e,0} = 0 + \end{cases} \\ + P_{\mathrm{ref},0} &= \text{clamp}(P_{e,0}, P^{\min}, P^{\max}) +\end{aligned} +``` + +When optional `qext`, `pfaref`, or `pref` signals are connected, initialization writes the recovered steady-state references to those signals. + +If $V_{\mathrm{ref0}}$ is omitted, set $V_{\mathrm{ref0}} = V_{T,0}$. + +```math +\begin{aligned} + P_{\mathrm{meas},0} &= P_{e,0} +\end{aligned} +``` + +Then evaluate the upstream algebraic chain: + +```math +\begin{aligned} + s_{\mathrm{dip},0} &= \text{outside}(V_{T,0}, V_{\mathrm{dip}}, V_{\mathrm{up}}) \\ + V_{\mathrm{err},0} &= \text{deadband2}(V_{\mathrm{ref0}} - V_{\mathrm{meas},0}, D_{\mathrm{bd1}}, D_{\mathrm{bd2}}) \\ + I_{\mathrm{qv},0} &= \text{clamp}(K_{\mathrm{qv}} V_{\mathrm{err},0}, I_{\mathrm{qinj}}^{\min}, I_{\mathrm{qinj}}^{\max}) \\ + Q_{\mathrm{ref},0} &= s_{\mathrm{pf}} P_{\mathrm{meas},0}\tan(\phi_{\mathrm{pf},0}^{\mathrm{ref}}) + s_{\mathrm{pf}}^{\mathrm{off}} Q_{\mathrm{ext},0} \\ + e_{Q,0} &= \text{clamp}(Q_{\mathrm{ref},0}, Q^{\min}, Q^{\max}) - Q_{\mathrm{gen},0} \\ + Q_{V,0} &= \dfrac{Q_{\mathrm{ref},0}}{V_{\mathrm{meas},0}^{\mathrm{safe}}} \\ + P_{\mathrm{ord},0} &= P_{\mathrm{ref},0} +\end{aligned} +``` + +Initialize the reactive-power PI output so its residual and zero-derivative anti-windup condition hold: + +```math +\begin{aligned} + V_{\mathrm{PIQ},0} &= \text{clamp}(K_{\mathrm{qp}} e_{Q,0} + x_{\mathrm{PIQ},0}, V^{\min}, V^{\max}) \\ + e_{\mathrm{PIV},0} &= s_V V_{\mathrm{PIQ},0} + s_V^{\mathrm{off}}Q_{\mathrm{ref},0} - V_{\mathrm{meas},0} +\end{aligned} +``` + +For an unsaturated zero-derivative start, require $e_{Q,0}=0$ for $x_{\mathrm{PIQ}}$ and choose or verify $e_{\mathrm{PIV},0}=0$ for $x_{\mathrm{PIV}}$. Then $x_{\mathrm{PIQ},0}=V_{\mathrm{PIQ},0}-K_{\mathrm{qp}}e_{Q,0}$; when $s_V=1$, set $V_{\mathrm{PIQ},0}=V_{\mathrm{meas},0}$, and when $s_V=0$, the selected $Q_{\mathrm{ref},0}$ signal must equal $V_{\mathrm{meas},0}$. For a saturated PI start, initialize the pre-limit PI output beyond the active limit in the direction that freezes the anti-windup gate. + +Finish by evaluating the current-limit and current-command algebraic residuals in priority order. At the command steps, use the power-flow current targets before final limiting: + +```math +\begin{aligned} + I_{\mathrm{qbase},0}^{\star} &= \dfrac{Q_{\mathrm{gen},0}}{V_{\mathrm{meas},0}^{\mathrm{safe}}} \\ + I_{\mathrm{p},0}^{\star} &= \dfrac{P_{\mathrm{ord},0}}{V_{\mathrm{meas},0}^{\mathrm{safe}}} +\end{aligned} +``` + +These current targets are in the terminal-voltage reference frame. REGCA rotates +the resulting `ipcmd` and `iqcmd` signals into network-frame current injection. + +For $s_{PQ}=0$, use: + +```math +\begin{aligned} +I_{\mathrm{q},0}^{\mathrm{circ}} +\rightarrow I_{\mathrm{q},0}^{\max} +\rightarrow I_{\mathrm{qbase},0} +\rightarrow I_{\mathrm{q},0}^{\mathrm{raw}} +\rightarrow I_{\mathrm{q},0}^{\mathrm{cmd}} +\rightarrow I_{\mathrm{p},0}^{\mathrm{circ}} +\rightarrow I_{\mathrm{p},0}^{\max} +\rightarrow I_{\mathrm{p},0}^{\mathrm{cmd}} +\end{aligned} +``` + +For $s_{PQ}=1$, use: + +```math +\begin{aligned} +I_{\mathrm{p},0}^{\mathrm{circ}} +\rightarrow I_{\mathrm{p},0}^{\max} +\rightarrow I_{\mathrm{p},0}^{\mathrm{cmd}} +\rightarrow I_{\mathrm{q},0}^{\mathrm{circ}} +\rightarrow I_{\mathrm{q},0}^{\max} +\rightarrow I_{\mathrm{qbase},0} +\rightarrow I_{\mathrm{q},0}^{\mathrm{raw}} +\rightarrow I_{\mathrm{q},0}^{\mathrm{cmd}} +\end{aligned} +``` + +After $I_{\mathrm{q},0}^{\max}$ and $I_{\mathrm{qbase},0}$ are known, initialize the voltage PI state from its output residual; the unsaturated zero-derivative start also requires the $e_{\mathrm{PIV},0}=0$ condition above: + +```math +\begin{aligned} + x_{\mathrm{PIV},0} &= I_{\mathrm{qbase},0} - K_{\mathrm{vp}} e_{\mathrm{PIV},0} +\end{aligned} +``` + +The current-circle variables use the nonnegative branch of the squared algebraic residuals, so the radicands must be nonnegative. A standard steady-state initialization assumes $s_{\mathrm{dip},0}=0$. Starts inside voltage-dip or overvoltage logic are not uniquely determined by these closed-form equations. + +## Model Outputs + +Output | Units | Description | Note +----------------|--------|-------------------------------------|------ +`iqcmd` | [p.u.] | Reactive-current command output | Converter base +`ipcmd` | [p.u.] | Active-current command output | Converter base +`vmeas` | [p.u.] | Filtered terminal voltage | +`pmeas` | [p.u.] | Filtered electrical power | +`piq` | [p.u.] | Reactive-power PI controller state | +`piv` | [p.u.] | Voltage PI controller state | +`qv` | [p.u.] | Reactive-current command lag state | +`pord` | [p.u.] | Filtered active-power order | +`qref` | [p.u.] | Selected reactive-power reference | +`sdip` | [binary] | Voltage-dip/overvoltage freeze indicator | +`iqmax` | [p.u.] | Final reactive-current upper limit | Converter base +`ipmax` | [p.u.] | Final active-current upper limit | Converter base +`iqv` | [p.u.] | Reactive-current injection candidate | Converter base +`vqctrl` | [p.u.] | Reactive-power control PI output | +`iqbase` | [p.u.] | Base reactive-current command | Converter base diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.cpp b/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.cpp new file mode 100644 index 000000000..170efdb20 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.cpp @@ -0,0 +1,27 @@ +/** + * @file Reecb.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the REECB electrical-control model. + */ + +#include "ReecbImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Reecb::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Reecb..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Reecb; + template class Reecb; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.hpp b/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.hpp new file mode 100644 index 000000000..44fd001ef --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.hpp @@ -0,0 +1,188 @@ +/** + * @file Reecb.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the REECB electrical-control model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + template + class SignalNode; + + namespace Converter + { + /// Internal variables of a `Reecb`. + enum class ReecbInternalVariables : size_t + { + VMEAS, ///< Filtered terminal voltage + PMEAS, ///< Filtered electrical power + XPIQ, ///< Reactive-power PI state + XPIV, ///< Voltage PI state + QV, ///< Reactive-current command lag state + PORD, ///< Filtered active-power order + VT, ///< Terminal voltage magnitude + VMEASSAFE, ///< Safe filtered terminal voltage + SDIP, ///< Voltage-dip/overvoltage freeze indicator + VERR, ///< Deadbanded voltage error + IQV, ///< Reactive-current injection candidate + QREF, ///< Selected reactive-power reference + EQ, ///< Reactive-power control error + VPIQ, ///< Reactive-power PI output + EPIV, ///< Voltage-control PI error + FPORD, ///< Active-power order derivative before ramp limiting + RPORD, ///< Ramp-rate-limited active-power derivative + IQCIRC, ///< Reactive-current limit from current circle + IPCIRC, ///< Active-current limit from current circle + IQMAX, ///< Final reactive-current upper limit + IPMAX, ///< Final active-current upper limit + IQBASE, ///< Base reactive-current command + IQRAW, ///< Raw reactive-current command + IQCMD, ///< Reactive-current command output + IPCMD, ///< Active-current command output + MAXIMUM, + }; + + /// External variables of a `Reecb`. + enum class ReecbExternalVariables : size_t + { + PE, ///< Electrical active-power signal + QGEN, ///< Reactive-power signal + QEXT, ///< External reactive-power command + PFAREF, ///< Power-factor angle reference + PREF, ///< External active-power reference + MAXIMUM, + }; + + template + class Reecb : 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::va_system_base_; + 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 = ReecbData; + using MonitorT = Model::VariableMonitor; + + Reecb(bus_type* bus); + Reecb(bus_type* bus, const model_data_type& data); + ~Reecb(); + + 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 initializeMonitor(); + void setDerivedParameters(); + + ScalarT toComponentBase(ScalarT value) const; + + ScalarT& Vr(); + ScalarT& Vi(); + + bus_type* bus_{nullptr}; + + RealT mva_base_{0}; + RealT PfFlag_{0}; + RealT VFlag_{0}; + RealT QFlag_{0}; + RealT Pqflag_{0}; + RealT Trv_{0}; + RealT Tp_{0}; + RealT Vref0_{0}; + RealT Vdip_{0}; + RealT Vup_{0}; + RealT dbd1_{0}; + RealT dbd2_{0}; + RealT kqv_{0}; + RealT Iql1_{0}; + RealT Iqh1_{0}; + RealT Qmax_{0}; + RealT Qmin_{0}; + RealT Kqp_{0}; + RealT Kqi_{0}; + RealT Vmax_{0}; + RealT Vmin_{0}; + RealT Kvp_{0}; + RealT Kvi_{0}; + RealT Tiq_{0}; + RealT Tpord_{0}; + RealT dPmax_{0}; + RealT dPmin_{0}; + RealT Pmax_{0}; + RealT Pmin_{0}; + RealT Imax_{0}; + RealT va_converter_base_{0}; + RealT pf_off_{1}; + RealT v_off_{1}; + RealT q_off_{1}; + + bool Vref0_given_{false}; + IdxT parameter_error_count_{0}; + + ScalarT qext_set_{0}; + ScalarT pe_set_{0}; + ScalarT qgen_set_{0}; + ScalarT pfaref_set_{0}; + ScalarT pref_set_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbData.hpp b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbData.hpp new file mode 100644 index 000000000..c421d33af --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbData.hpp @@ -0,0 +1,100 @@ +/** + * @file ReecbData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the REECB electrical-control model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + /// Parameter keys for the REECB electrical-control model. + enum class ReecbParameters + { + mva, ///< Model MVA base + PfFlag, ///< Power-factor control flag + VFlag, ///< Voltage-control mode flag + QFlag, ///< Reactive-power control flag + Pqflag, ///< P/Q current-priority flag + Trv, ///< Voltage-measurement filter time constant + Tp, ///< Electrical-power measurement filter time constant + Vref0, ///< Outer-loop voltage reference + Vdip, ///< Low-voltage dip threshold + Vup, ///< High-voltage threshold + dbd1, ///< Overvoltage deadband threshold + dbd2, ///< Undervoltage deadband threshold + kqv, ///< Reactive-current injection gain + Iql1, ///< Minimum reactive-current injection limit + Iqh1, ///< Maximum reactive-current injection limit + Qmax, ///< Maximum reactive-power control limit + Qmin, ///< Minimum reactive-power control limit + Kqp, ///< Reactive-power proportional gain + Kqi, ///< Reactive-power integral gain + Vmax, ///< Maximum voltage-control limit + Vmin, ///< Minimum voltage-control limit + Kvp, ///< Voltage-control proportional gain + Kvi, ///< Voltage-control integral gain + Tiq, ///< Reactive-current command lag time constant + Tpord, ///< Active-power order filter time constant + dPmax, ///< Positive active-power ramp-rate limit + dPmin, ///< Negative active-power ramp-rate limit + Pmax, ///< Maximum active-power order limit + Pmin, ///< Minimum active-power order limit + Imax ///< Maximum converter current + }; + + /// Ports for the REECB electrical-control model. + enum class ReecbPorts + { + bus, ///< Terminal bus ID + pe, ///< Electrical active-power signal ID + qgen, ///< Reactive-power signal ID + qext, ///< Optional reactive-power command signal ID + pfaref, ///< Optional power-factor angle reference signal ID + pref, ///< Optional active-power reference signal ID + iqcmd, ///< Reactive-current command output signal ID + ipcmd ///< Active-current command output signal ID + }; + + /// Variables available through the monitor interface. + enum class ReecbMonitorableVariables + { + iqcmd, ///< Reactive-current command output + ipcmd, ///< Active-current command output + vmeas, ///< Filtered terminal voltage + pmeas, ///< Filtered electrical power + piq, ///< Reactive-power PI controller state + piv, ///< Voltage PI controller state + qv, ///< Reactive-current command lag state + pord, ///< Filtered active-power order + qref, ///< Selected reactive-power reference + sdip, ///< Voltage-dip/overvoltage freeze indicator + iqmax, ///< Final reactive-current upper limit + ipmax, ///< Final active-current upper limit + iqv, ///< Reactive-current injection candidate + vqctrl, ///< Reactive-power control PI output + iqbase ///< Base reactive-current command + }; + + template + struct ReecbData : public ComponentData + { + ReecbData() = default; + + using Parameters = ReecbParameters; + using Ports = ReecbPorts; + using MonitorableVariables = ReecbMonitorableVariables; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbDependencyTracking.cpp new file mode 100644 index 000000000..4ef115e1d --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file ReecbDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the REECB electrical-control model. + */ + +#include "ReecbImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Reecb::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Reecb..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Reecb; + template class Reecb; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbEnzyme.cpp new file mode 100644 index 000000000..5c077b3ad --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbEnzyme.cpp @@ -0,0 +1,92 @@ +/** + * @file ReecbEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the REECB electrical-control model. + */ + +#include + +#include "ReecbImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Reecb::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Reecb..." << 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::Converter::Reecb; + 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(), + wb_.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 Reecb; + template class Reecb; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbImpl.hpp new file mode 100644 index 000000000..d484971ec --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbImpl.hpp @@ -0,0 +1,721 @@ +/** + * @file ReecbImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the REECB electrical-control model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + using Log = ::GridKit::Utilities::Logger; + + template + Reecb::Reecb(bus_type* bus) + : bus_(bus) + { + size_ = static_cast(ReecbInternalVariables::MAXIMUM); + } + + template + Reecb::Reecb(bus_type* bus, const model_data_type& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initModelParams(data); + initializeMonitor(); + size_ = static_cast(ReecbInternalVariables::MAXIMUM); + } + + template + Reecb::~Reecb() + { + } + + template + ScalarT& Reecb::Vr() + { + return bus_->Vr(); + } + + template + ScalarT& Reecb::Vi() + { + return bus_->Vi(); + } + + template + void Reecb::setDerivedParameters() + { + va_converter_base_ = mva_base_ * static_cast(1.0e6); + pf_off_ = ONE - PfFlag_; + v_off_ = ONE - VFlag_; + q_off_ = ONE - QFlag_; + } + + template + ScalarT Reecb::toComponentBase(ScalarT value) const + { + return value * va_system_base_ / va_converter_base_; + } + + template + void Reecb::initModelParams(const model_data_type& data) + { + using Params = typename model_data_type::Parameters; + + parameter_error_count_ = 0; + Vref0_given_ = false; + + auto loadRequiredReal = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Reecb: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* real_value = std::get_if(&value)) + { + target = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + target = static_cast(*index_value); + } + else + { + Log::error() << "Reecb: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + auto loadRequiredSwitch = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Reecb: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* bool_value = std::get_if(&value)) + { + target = ZERO; + if (*bool_value) + { + target = ONE; + } + } + else if (const auto* index_value = std::get_if(&value); + index_value && (*index_value == 0 || *index_value == 1)) + { + target = static_cast(*index_value); + } + else if (const auto* real_value = std::get_if(&value); + real_value && (*real_value == ZERO || *real_value == ONE) ) + { + target = *real_value; + } + else + { + Log::error() << "Reecb: parameter '" << name << "' must be bool or 0/1\n"; + ++parameter_error_count_; + } + }; + + loadRequiredReal(Params::mva, mva_base_, "mva"); + loadRequiredSwitch(Params::PfFlag, PfFlag_, "PfFlag"); + loadRequiredSwitch(Params::VFlag, VFlag_, "VFlag"); + loadRequiredSwitch(Params::QFlag, QFlag_, "QFlag"); + loadRequiredSwitch(Params::Pqflag, Pqflag_, "Pqflag"); + loadRequiredReal(Params::Trv, Trv_, "Trv"); + loadRequiredReal(Params::Tp, Tp_, "Tp"); + if (data.parameters.contains(Params::Vref0)) + { + loadRequiredReal(Params::Vref0, Vref0_, "Vref0"); + Vref0_given_ = true; + } + loadRequiredReal(Params::Vdip, Vdip_, "Vdip"); + loadRequiredReal(Params::Vup, Vup_, "Vup"); + loadRequiredReal(Params::dbd1, dbd1_, "dbd1"); + loadRequiredReal(Params::dbd2, dbd2_, "dbd2"); + loadRequiredReal(Params::kqv, kqv_, "kqv"); + loadRequiredReal(Params::Iql1, Iql1_, "Iql1"); + loadRequiredReal(Params::Iqh1, Iqh1_, "Iqh1"); + loadRequiredReal(Params::Qmax, Qmax_, "Qmax"); + loadRequiredReal(Params::Qmin, Qmin_, "Qmin"); + loadRequiredReal(Params::Kqp, Kqp_, "Kqp"); + loadRequiredReal(Params::Kqi, Kqi_, "Kqi"); + loadRequiredReal(Params::Vmax, Vmax_, "Vmax"); + loadRequiredReal(Params::Vmin, Vmin_, "Vmin"); + loadRequiredReal(Params::Kvp, Kvp_, "Kvp"); + loadRequiredReal(Params::Kvi, Kvi_, "Kvi"); + loadRequiredReal(Params::Tiq, Tiq_, "Tiq"); + loadRequiredReal(Params::Tpord, Tpord_, "Tpord"); + loadRequiredReal(Params::dPmax, dPmax_, "dPmax"); + loadRequiredReal(Params::dPmin, dPmin_, "dPmin"); + loadRequiredReal(Params::Pmax, Pmax_, "Pmax"); + loadRequiredReal(Params::Pmin, Pmin_, "Pmin"); + loadRequiredReal(Params::Imax, Imax_, "Imax"); + setDerivedParameters(); + } + + template + const Model::VariableMonitorBase* Reecb::getMonitor() const + { + return monitor_.get(); + } + + template + void Reecb::initializeMonitor() + { + using Variable = typename model_data_type::MonitorableVariables; + auto index = [](ReecbInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::iqcmd, [this, index] + { return y_[index(ReecbInternalVariables::IQCMD)]; }); + monitor_->set(Variable::ipcmd, [this, index] + { return y_[index(ReecbInternalVariables::IPCMD)]; }); + monitor_->set(Variable::vmeas, [this, index] + { return y_[index(ReecbInternalVariables::VMEAS)]; }); + monitor_->set(Variable::pmeas, [this, index] + { return y_[index(ReecbInternalVariables::PMEAS)]; }); + monitor_->set(Variable::piq, [this, index] + { return y_[index(ReecbInternalVariables::XPIQ)]; }); + monitor_->set(Variable::piv, [this, index] + { return y_[index(ReecbInternalVariables::XPIV)]; }); + monitor_->set(Variable::qv, [this, index] + { return y_[index(ReecbInternalVariables::QV)]; }); + monitor_->set(Variable::pord, [this, index] + { return y_[index(ReecbInternalVariables::PORD)]; }); + monitor_->set(Variable::qref, [this, index] + { return y_[index(ReecbInternalVariables::QREF)]; }); + monitor_->set(Variable::sdip, [this, index] + { return y_[index(ReecbInternalVariables::SDIP)]; }); + monitor_->set(Variable::iqmax, [this, index] + { return y_[index(ReecbInternalVariables::IQMAX)]; }); + monitor_->set(Variable::ipmax, [this, index] + { return y_[index(ReecbInternalVariables::IPMAX)]; }); + monitor_->set(Variable::iqv, [this, index] + { return y_[index(ReecbInternalVariables::IQV)]; }); + monitor_->set(Variable::vqctrl, [this, index] + { return y_[index(ReecbInternalVariables::VPIQ)]; }); + monitor_->set(Variable::iqbase, [this, index] + { return y_[index(ReecbInternalVariables::IQBASE)]; }); + } + + template + int Reecb::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Reecb::allocate() + { + size_ = static_cast(ReecbInternalVariables::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(ReecbExternalVariables::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(ReecbInternalVariables::IQCMD)], + &(this->getVariableIndex(static_cast(ReecbInternalVariables::IQCMD)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(ReecbInternalVariables::IPCMD)], + &(this->getVariableIndex(static_cast(ReecbInternalVariables::IPCMD)))); + } + + return 0; + } + + template + int Reecb::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Reecb: " << message << '\n'; + ret += 1; + } + }; + + if (bus_ == nullptr) + { + Log::error() << "Reecb: bus pointer is null\n"; + ret += 1; + } + + check(mva_base_ > ZERO, "mva must be positive"); + check(va_converter_base_ > ZERO, "converter VA base must be positive"); + check(PfFlag_ == ZERO || PfFlag_ == ONE, "PfFlag must be 0 or 1"); + check(VFlag_ == ZERO || VFlag_ == ONE, "VFlag must be 0 or 1"); + check(QFlag_ == ZERO || QFlag_ == ONE, "QFlag must be 0 or 1"); + check(Pqflag_ == ZERO || Pqflag_ == ONE, "Pqflag must be 0 or 1"); + check(Trv_ >= ZERO, "Trv must be non-negative"); + check(Tp_ >= ZERO, "Tp must be non-negative"); + check(Vdip_ < Vup_, "Vdip must be less than Vup"); + check(dbd1_ <= ZERO && ZERO <= dbd2_, "dbd1 <= 0 <= dbd2 is required"); + check(Iql1_ <= Iqh1_, "Iql1 must be less than or equal to Iqh1"); + check(Qmin_ <= Qmax_, "Qmin must be less than or equal to Qmax"); + check(Vmin_ <= Vmax_, "Vmin must be less than or equal to Vmax"); + check(Tiq_ > ZERO, "Tiq must be positive"); + check(Tpord_ > ZERO, "Tpord must be positive"); + check(dPmin_ < ZERO && ZERO < dPmax_, "dPmin < 0 < dPmax is required"); + check(Pmin_ <= Pmax_, "Pmin must be less than or equal to Pmax"); + check(Imax_ >= ZERO, "Imax must be non-negative"); + + const bool has_pe = signals_.template isAttached(); + const bool has_qgen = signals_.template isAttached(); + + if (has_pe != has_qgen) + { + Log::error() << "Reecb: pe and qgen must be connected together\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Reecb: pe signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Reecb: qgen signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Reecb: qext signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Reecb: pfaref signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Reecb: pref signal attached with no linked source\n"; + ret += 1; + } + + return ret; + } + + template + int Reecb::initialize() + { + if (parameter_error_count_ > 0 || verify() > 0) + { + Log::error() << "Reecb: cannot initialize with invalid configuration\n"; + return 1; + } + + const auto VMEAS = static_cast(ReecbInternalVariables::VMEAS); + const auto PMEAS = static_cast(ReecbInternalVariables::PMEAS); + const auto XPIQ = static_cast(ReecbInternalVariables::XPIQ); + const auto XPIV = static_cast(ReecbInternalVariables::XPIV); + const auto QV = static_cast(ReecbInternalVariables::QV); + const auto PORD = static_cast(ReecbInternalVariables::PORD); + const auto VT = static_cast(ReecbInternalVariables::VT); + const auto VMEASSAFE = static_cast(ReecbInternalVariables::VMEASSAFE); + const auto SDIP = static_cast(ReecbInternalVariables::SDIP); + const auto VERR = static_cast(ReecbInternalVariables::VERR); + const auto IQV = static_cast(ReecbInternalVariables::IQV); + const auto QREF = static_cast(ReecbInternalVariables::QREF); + const auto EQ = static_cast(ReecbInternalVariables::EQ); + const auto VPIQ = static_cast(ReecbInternalVariables::VPIQ); + const auto EPIV = static_cast(ReecbInternalVariables::EPIV); + const auto FPORD = static_cast(ReecbInternalVariables::FPORD); + const auto RPORD = static_cast(ReecbInternalVariables::RPORD); + const auto IQCIRC = static_cast(ReecbInternalVariables::IQCIRC); + const auto IPCIRC = static_cast(ReecbInternalVariables::IPCIRC); + const auto IQMAX = static_cast(ReecbInternalVariables::IQMAX); + const auto IPMAX = static_cast(ReecbInternalVariables::IPMAX); + const auto IQBASE = static_cast(ReecbInternalVariables::IQBASE); + const auto IQRAW = static_cast(ReecbInternalVariables::IQRAW); + const auto IQCMD = static_cast(ReecbInternalVariables::IQCMD); + const auto IPCMD = static_cast(ReecbInternalVariables::IPCMD); + + const ScalarT vr = Vr(); + const ScalarT vi = Vi(); + + y_[VT] = std::sqrt(vr * vr + vi * vi); + if (!Vref0_given_) + { + Vref0_ = static_cast(y_[VT]); + } + + y_[VMEAS] = y_[VT]; + y_[VMEASSAFE] = Math::max(y_[VMEAS], static_cast(0.01)); + + const ScalarT ipcmd0 = y_[IPCMD]; + const ScalarT iqcmd0 = y_[IQCMD]; + const ScalarT pe0 = ipcmd0 * y_[VMEASSAFE]; + const ScalarT qgen0 = iqcmd0 * y_[VMEASSAFE]; + + const ScalarT qref0 = qgen0; + const ScalarT pref0 = Math::clamp(pe0, Pmin_, Pmax_); + ScalarT pfaref0{ZERO}; + if (std::abs(static_cast(pe0)) > ZERO) + { + pfaref0 = static_cast(std::atan(static_cast(qref0 / pe0))); + } + + pe_set_ = pe0; + qgen_set_ = qgen0; + qext_set_ = qref0; + pfaref_set_ = pfaref0; + pref_set_ = pref0; + + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(qref0); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(pfaref0); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(pref0); + } + + y_[PMEAS] = pe0; + const ScalarT qref_pf = PfFlag_ * y_[PMEAS] * std::tan(pfaref_set_); + const ScalarT qref_ext = pf_off_ * qext_set_; + + y_[SDIP] = Math::outside(y_[VT], Vdip_, Vup_); + y_[VERR] = Math::deadband2(Vref0_ - y_[VMEAS], dbd1_, dbd2_); + y_[IQV] = Math::clamp(kqv_ * y_[VERR], Iql1_, Iqh1_); + y_[QREF] = qref_pf + qref_ext; + y_[EQ] = Math::clamp(y_[QREF], Qmin_, Qmax_) - qgen_set_; + y_[QV] = y_[QREF] / y_[VMEASSAFE]; + y_[PORD] = pref_set_; + y_[FPORD] = ZERO; + y_[RPORD] = ZERO; + + static constexpr RealT kSmallTol = static_cast(1.0e-9); + static constexpr RealT kSatMargin = static_cast(0.1); + + auto steadyPiOutput = [&](const ScalarT interior_target, + const ScalarT error, + const ScalarT integral_rate, + const ScalarT lower, + const ScalarT upper) -> ScalarT + { + const RealT error_value = static_cast(error); + const RealT rate_value = static_cast(integral_rate); + + if (std::abs(error_value) <= kSmallTol || std::abs(rate_value) <= kSmallTol) + { + return interior_target; + } + + if (rate_value > ZERO) + { + return upper + static_cast(kSatMargin); + } + + return lower - static_cast(kSatMargin); + }; + + const ScalarT xpiq_rate = Kqi_ * y_[EQ]; + const ScalarT vpiq_target = VFlag_ * y_[VMEAS] + v_off_ * y_[QREF]; + const ScalarT vpiq_arg = steadyPiOutput(vpiq_target, + y_[EQ], + xpiq_rate, + static_cast(Vmin_), + static_cast(Vmax_)); + y_[VPIQ] = Math::clamp(vpiq_arg, Vmin_, Vmax_); + y_[EPIV] = VFlag_ * y_[VPIQ] + v_off_ * y_[QREF] - y_[VMEAS]; + y_[XPIQ] = vpiq_arg - Kqp_ * y_[EQ]; + + const ScalarT iqbase_target = qgen0 / y_[VMEASSAFE]; + const ScalarT xpiv_rate = Kvi_ * y_[EPIV]; + const ScalarT ip_star = y_[PORD] / y_[VMEASSAFE]; + + auto initializeReactiveBase = [&]() + { + const ScalarT piv_lower = -y_[IQMAX]; + const ScalarT piv_upper = y_[IQMAX]; + const ScalarT piv_arg = steadyPiOutput(iqbase_target, + y_[EPIV], + xpiv_rate, + piv_lower, + piv_upper); + y_[IQBASE] = Math::clamp(piv_arg, -y_[IQMAX], y_[IQMAX]); + y_[XPIV] = piv_arg - Kvp_ * y_[EPIV]; + }; + + auto initializeReactiveCommand = [&]() + { + y_[IQRAW] = QFlag_ * y_[IQBASE] + q_off_ * y_[QV] + y_[SDIP] * y_[IQV]; + y_[IQCMD] = Math::clamp(y_[IQRAW], -y_[IQMAX], y_[IQMAX]); + }; + + if (Pqflag_ == ZERO) + { + y_[IQCIRC] = Imax_; + y_[IQMAX] = Imax_; + initializeReactiveBase(); + initializeReactiveCommand(); + + const ScalarT ip_radicand = Imax_ * Imax_ - y_[IQCMD] * y_[IQCMD]; + if (static_cast(ip_radicand) < ZERO) + { + Log::error() << "Reecb: initial active-current circle radicand is negative\n"; + return 1; + } + y_[IPCIRC] = std::sqrt(ip_radicand); + y_[IPMAX] = y_[IPCIRC]; + y_[IPCMD] = Math::clamp(ip_star, ZERO, y_[IPMAX]); + } + else + { + y_[IPCIRC] = Imax_; + y_[IPMAX] = Imax_; + y_[IPCMD] = Math::clamp(ip_star, ZERO, y_[IPMAX]); + + const ScalarT iq_radicand = Imax_ * Imax_ - y_[IPCMD] * y_[IPCMD]; + if (static_cast(iq_radicand) < ZERO) + { + Log::error() << "Reecb: initial reactive-current circle radicand is negative\n"; + return 1; + } + y_[IQCIRC] = std::sqrt(iq_radicand); + y_[IQMAX] = y_[IQCIRC]; + initializeReactiveBase(); + initializeReactiveCommand(); + } + + std::fill(yp_.begin(), yp_.end(), ZERO); + + return 0; + } + + template + int Reecb::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(ReecbInternalVariables::VMEAS)] = (Trv_ > ZERO); + tag_[static_cast(ReecbInternalVariables::PMEAS)] = (Tp_ > ZERO); + tag_[static_cast(ReecbInternalVariables::XPIQ)] = true; + tag_[static_cast(ReecbInternalVariables::XPIV)] = true; + tag_[static_cast(ReecbInternalVariables::QV)] = true; + tag_[static_cast(ReecbInternalVariables::PORD)] = true; + return 0; + } + + template + __attribute__((always_inline)) inline int Reecb::evaluateInternalResidual( + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + ScalarT* f) + { + const auto VMEAS = static_cast(ReecbInternalVariables::VMEAS); + const auto PMEAS = static_cast(ReecbInternalVariables::PMEAS); + const auto XPIQ = static_cast(ReecbInternalVariables::XPIQ); + const auto XPIV = static_cast(ReecbInternalVariables::XPIV); + const auto QV = static_cast(ReecbInternalVariables::QV); + const auto PORD = static_cast(ReecbInternalVariables::PORD); + const auto VT = static_cast(ReecbInternalVariables::VT); + const auto VMEASSAFE = static_cast(ReecbInternalVariables::VMEASSAFE); + const auto SDIP = static_cast(ReecbInternalVariables::SDIP); + const auto VERR = static_cast(ReecbInternalVariables::VERR); + const auto IQV = static_cast(ReecbInternalVariables::IQV); + const auto QREF = static_cast(ReecbInternalVariables::QREF); + const auto EQ = static_cast(ReecbInternalVariables::EQ); + const auto VPIQ = static_cast(ReecbInternalVariables::VPIQ); + const auto EPIV = static_cast(ReecbInternalVariables::EPIV); + const auto FPORD = static_cast(ReecbInternalVariables::FPORD); + const auto RPORD = static_cast(ReecbInternalVariables::RPORD); + const auto IQCIRC = static_cast(ReecbInternalVariables::IQCIRC); + const auto IPCIRC = static_cast(ReecbInternalVariables::IPCIRC); + const auto IQMAX = static_cast(ReecbInternalVariables::IQMAX); + const auto IPMAX = static_cast(ReecbInternalVariables::IPMAX); + const auto IQBASE = static_cast(ReecbInternalVariables::IQBASE); + const auto IQRAW = static_cast(ReecbInternalVariables::IQRAW); + const auto IQCMD = static_cast(ReecbInternalVariables::IQCMD); + const auto IPCMD = static_cast(ReecbInternalVariables::IPCMD); + + const auto PE = static_cast(ReecbExternalVariables::PE); + const auto QGEN = static_cast(ReecbExternalVariables::QGEN); + const auto QEXT = static_cast(ReecbExternalVariables::QEXT); + const auto PFAREF = static_cast(ReecbExternalVariables::PFAREF); + const auto PREF = static_cast(ReecbExternalVariables::PREF); + + const ScalarT vr = wb[0]; + const ScalarT vi = wb[1]; + + const ScalarT pe = ws[PE]; + const ScalarT qgen = ws[QGEN]; + const ScalarT qext = ws[QEXT]; + const ScalarT pfaref = ws[PFAREF]; + const ScalarT pref = ws[PREF]; + + const ScalarT not_sdip = ONE - y[SDIP]; + const ScalarT xpiq_arg = Kqp_ * y[EQ] + y[XPIQ]; + const ScalarT xpiv_arg = Kvp_ * y[EPIV] + y[XPIV]; + const ScalarT xpiq_limited_rate = Math::antiwindup(xpiq_arg, Kqi_ * y[EQ], Vmin_, Vmax_); + const ScalarT xpiv_limited_rate = Math::antiwindup(xpiv_arg, Kvi_ * y[EPIV], -y[IQMAX], y[IQMAX]); + const ScalarT pord_limited_rate = Math::antiwindup(y[PORD], y[RPORD], Pmin_, Pmax_); + const ScalarT qref_pf = PfFlag_ * y[PMEAS] * std::tan(pfaref); + const ScalarT qref_ext = pf_off_ * qext; + const ScalarT epiv_target = VFlag_ * y[VPIQ] + v_off_ * y[QREF]; + const ScalarT iqcirc_current_sq = Pqflag_ * y[IPCMD] * y[IPCMD]; + const ScalarT ipcirc_current_sq = (ONE - Pqflag_) * y[IQCMD] * y[IQCMD]; + const ScalarT iqmax_target = (ONE - Pqflag_) * Imax_ + Pqflag_ * y[IQCIRC]; + const ScalarT ipmax_target = Pqflag_ * Imax_ + (ONE - Pqflag_) * y[IPCIRC]; + const ScalarT iqraw_target = QFlag_ * y[IQBASE] + q_off_ * y[QV] + y[SDIP] * y[IQV]; + + f[VMEAS] = -Trv_ * yp[VMEAS] - y[VMEAS] + y[VT]; + f[PMEAS] = -Tp_ * yp[PMEAS] - y[PMEAS] + pe; + f[XPIQ] = -yp[XPIQ] + not_sdip * xpiq_limited_rate; + f[XPIV] = -yp[XPIV] + not_sdip * xpiv_limited_rate; + f[QV] = -Tiq_ * yp[QV] - not_sdip * y[QV] + not_sdip * y[QREF] / y[VMEASSAFE]; + f[PORD] = -yp[PORD] + not_sdip * pord_limited_rate; + + f[VT] = -y[VT] * y[VT] + vr * vr + vi * vi; + f[VMEASSAFE] = -y[VMEASSAFE] + Math::max(y[VMEAS], static_cast(0.01)); + f[SDIP] = -y[SDIP] + Math::outside(y[VT], Vdip_, Vup_); + f[VERR] = -y[VERR] + Math::deadband2(Vref0_ - y[VMEAS], dbd1_, dbd2_); + f[IQV] = -y[IQV] + Math::clamp(kqv_ * y[VERR], Iql1_, Iqh1_); + f[QREF] = -y[QREF] + qref_pf + qref_ext; + f[EQ] = -y[EQ] + Math::clamp(y[QREF], Qmin_, Qmax_) - qgen; + f[VPIQ] = -y[VPIQ] + Math::clamp(xpiq_arg, Vmin_, Vmax_); + f[EPIV] = -y[EPIV] + epiv_target - y[VMEAS]; + f[FPORD] = -Tpord_ * y[FPORD] + pref - y[PORD]; + f[RPORD] = -y[RPORD] + Math::clamp(y[FPORD], dPmin_, dPmax_); + + f[IQCIRC] = -y[IQCIRC] * y[IQCIRC] + Imax_ * Imax_ - iqcirc_current_sq; + f[IPCIRC] = -y[IPCIRC] * y[IPCIRC] + Imax_ * Imax_ - ipcirc_current_sq; + f[IQMAX] = -y[IQMAX] + iqmax_target; + f[IPMAX] = -y[IPMAX] + ipmax_target; + f[IQBASE] = -y[IQBASE] + Math::clamp(xpiv_arg, -y[IQMAX], y[IQMAX]); + f[IQRAW] = -y[IQRAW] + iqraw_target; + f[IQCMD] = -y[IQCMD] + Math::clamp(y[IQRAW], -y[IQMAX], y[IQMAX]); + f[IPCMD] = -y[IPCMD] + Math::clamp(y[PORD] / y[VMEASSAFE], ZERO, y[IPMAX]); + + return 0; + } + + template + int Reecb::evaluateResidual() + { + const auto PE = static_cast(ReecbExternalVariables::PE); + const auto QGEN = static_cast(ReecbExternalVariables::QGEN); + const auto QEXT = static_cast(ReecbExternalVariables::QEXT); + const auto PFAREF = static_cast(ReecbExternalVariables::PFAREF); + const auto PREF = static_cast(ReecbExternalVariables::PREF); + + ws_[PE] = pe_set_; + ws_[QGEN] = qgen_set_; + ws_[QEXT] = qext_set_; + ws_[PFAREF] = pfaref_set_; + ws_[PREF] = pref_set_; + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + if (signals_.template isAttached()) + { + ws_[PE] = toComponentBase( + signals_.template readExternalVariable()); + ws_indices_[PE] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[QGEN] = toComponentBase( + signals_.template readExternalVariable()); + ws_indices_[QGEN] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[QEXT] = signals_.template readExternalVariable(); + ws_indices_[QEXT] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[PFAREF] = signals_.template readExternalVariable(); + ws_indices_[PFAREF] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[PREF] = signals_.template readExternalVariable(); + ws_indices_[PREF] = + signals_.template readExternalVariableIndex(); + } + + wb_[0] = Vr(); + wb_[1] = Vi(); + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + return 0; + } + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index b33849140..4918b47f4 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -145,6 +145,7 @@ are specified: `Gensal` | 5th order salient-pole machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed`, `Eqp`, `psidp`, `psiqpp`, `psidpp`, `vd`, `vq`, `te`, `id`, `iq` `GenClassical` | the classical machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Xdp`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega` `Regca` | WECC REGCA renewable generator/converter model | `bus`, `ipcmd`\*, `iqcmd`\*, `ibranchr`\*, `ibranchi`\*, `pbranch`\*, `qbranch`\* | `P0`, `Q0`, `mva`, `Tg`, `TM`, `Rqmax`, `Rqmin`, `Rpmax`, `sL`, `IL1`, `VL0`, `VL1`, `VA0`, `VA1`, `Vhvmax` | `ir`, `ii`, `p`, `q`, `vt`, `vm`, `ip`, `iq`, `iqextra`, `il`, `lp`, `up` + `Reecb` | the REECB renewable electrical-control model | `bus`, `pe`\*, `qgen`\*, `qext`\*, `pfaref`\*, `pref`\*, `iqcmd`, `ipcmd` | `mva`, `PfFlag`, `VFlag`, `QFlag`, `Pqflag`, `Trv`, `Tp`, `Vref0`, `Vdip`, `Vup`, `dbd1`, `dbd2`, `kqv`, `Iql1`, `Iqh1`, `Qmax`, `Qmin`, `Kqp`, `Kqi`, `Vmax`, `Vmin`, `Kvp`, `Kvi`, `Tiq`, `Tpord`, `dPmax`, `dPmin`, `Pmax`, `Pmin`, `Imax` | `iqcmd`, `ipcmd`, `vmeas`, `pmeas`, `piq`, `piv`, `qv`, `pord`, `qref`, `sdip`, `iqmax`, `ipmax`, `iqv`, `vqctrl`, `iqbase` `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` `SexsPti` | the SEXS-PTI simplified exciter model | `bus`, `efd`, `vs`\* | `Ta`, `Tb`, `Te`, `K`, `Efdmax`, `Efdmin` | `efd` @@ -153,7 +154,8 @@ are specified: `BusToSignalAdapter` | signal adapter component for a bus | `bus`, `vr`, `vi`, `ir`, `ii` | | Ports marked with \* are optional and, if missing, will be assumed to be -connected to a constant value. This list is subject to change. +connected to a constant value. +This list is subject to change. For `Branch`, `tap` and `phase` are optional parameters. If omitted, `tap` defaults to `1.0` and `phase` defaults to `0.0` radians. Bus `bus1` is the tap diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 00b0e73ab..9098a950a 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -152,7 +152,7 @@ namespace GridKit addComponent(adapter); } - // Add REGCA converters + // Add REGCA converters before REECB so REGCA initializes the current commands first. for (const auto& regcadata : data.regca) { using DataT = typename SystemModelData::RegcaDataT; @@ -210,6 +210,65 @@ namespace GridKit addComponent(regca); } + // Add REECB electrical controllers + for (const auto& reecbdata : data.reecb) + { + using DataT = typename SystemModelData::ReecbDataT; + + IdxT bus_index = 0; + if (reecbdata.ports.contains(DataT::Ports::bus)) + { + bus_index = reecbdata.ports.at(DataT::Ports::bus); + } + + auto* reecb = new Reecb(getBus(bus_index), reecbdata); + + if (reecbdata.ports.contains(DataT::Ports::pe)) + { + const IdxT signal = reecbdata.ports.at(DataT::Ports::pe); + reecb->getSignals().template attachSignalNode( + getSignal(signal)); + } + if (reecbdata.ports.contains(DataT::Ports::qgen)) + { + const IdxT signal = reecbdata.ports.at(DataT::Ports::qgen); + reecb->getSignals().template attachSignalNode( + getSignal(signal)); + } + if (reecbdata.ports.contains(DataT::Ports::qext)) + { + const IdxT signal = reecbdata.ports.at(DataT::Ports::qext); + reecb->getSignals().template attachSignalNode( + getSignal(signal)); + } + if (reecbdata.ports.contains(DataT::Ports::pfaref)) + { + const IdxT signal = reecbdata.ports.at(DataT::Ports::pfaref); + reecb->getSignals().template attachSignalNode( + getSignal(signal)); + } + if (reecbdata.ports.contains(DataT::Ports::pref)) + { + const IdxT signal = reecbdata.ports.at(DataT::Ports::pref); + reecb->getSignals().template attachSignalNode( + getSignal(signal)); + } + if (reecbdata.ports.contains(DataT::Ports::iqcmd)) + { + const IdxT signal = reecbdata.ports.at(DataT::Ports::iqcmd); + reecb->getSignals().template assignSignalNode( + getSignal(signal)); + } + if (reecbdata.ports.contains(DataT::Ports::ipcmd)) + { + const IdxT signal = reecbdata.ports.at(DataT::Ports::ipcmd); + reecb->getSignals().template assignSignalNode( + getSignal(signal)); + } + + addComponent(reecb); + } + // Add branches for (const auto& branchdata : data.branch) { diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index 20501d169..59e5cbd28 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -40,6 +41,7 @@ namespace GridKit using BusToSignalAdapterDataT = BusToSignalAdapterData; using BusFaultDataT = BusFaultData; using RegcaDataT = Converter::RegcaData; + using ReecbDataT = Converter::ReecbData; using Tgov1DataT = Governor::Tgov1Data; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; @@ -94,6 +96,7 @@ namespace GridKit std::vector branch; ///< Branches within the model std::vector bus_fault; ///< Bus faults within the model std::vector regca; ///< REGCA converter instances within the model + std::vector reecb; ///< REECB electrical controllers within the model std::vector genrou; ///< GENROU instances within the model std::vector gensal; ///< GENSAL instances within the model std::vector genclassical; ///< Classical generator instances within the model diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index b215c0c86..bcc46e3ed 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -142,6 +142,12 @@ namespace GridKit raw_component.get_to(regca); sm.regca.push_back(regca); } + else if (kind == "Reecb") + { + typename SystemModelData::ReecbDataT reecb; + raw_component.get_to(reecb); + sm.reecb.push_back(reecb); + } else if (kind == "Tgov1") { typename SystemModelData::Tgov1DataT gov; diff --git a/docs/Figures/PhasorDynamics_REECB_Diagram.png b/docs/Figures/PhasorDynamics_REECB_Diagram.png new file mode 100644 index 0000000000000000000000000000000000000000..c3aaa1478e941c1b11a4b616843b840270ab4720 GIT binary patch literal 77117 zcmd43g;$hq)IK_NcS?hUA|fbK(k%mufTYqbAc%CQq=ck^NJ}cE0wN_PASo#g(jcL9 zojttYx6b*kbN+y{*030cndiA<-+N#C+Si_NHPyQ$gbaiT1cC&mB(H%$V96p7m^%2l z@Jc=Z?`8Of;iPd_7V)8nX$AhkwvxFogFuwW5S^Idz~2PVl=Pes2-0Ts4@SE~fjI(k zo`;f`(R4RlpLBae)q0G#)!QFqDEf11`KbRxpG5q?PPutwwd`rPMqNw-#_r^)LX_7j_ULCj&>h{Iu z!U^<}Om5Q8xYlAITx*$GV(@4qUrZ16Ud`|h95M1^MyO5bXsYjePxj~BEw|{R{kl8+ z_-0zSuM6CNfb3!zrftSI8%@ilOGldRgqlsusmB0`#8L;ld!pgBH;P55JVt zkDuIygCNHaz(fz8APAEke%;C$Su?2`heuYw^}(vYym~NKgujDGI#1ULFWf~bT%AoG zpB-tb%NN@`>)&A>$kSCez7dNgh(lu3v)i6*D&|mRbA2=Cm!Vc548cdKt@iSM0Q6Wu0XyNOrq!dyvYN{&3jf+*V9##cLK=Q1cOTNmQ;sWWX{Ak7$hjYeNY#} z)73f^nu{0IBrRu%8pxIogCX3z6xAl1GWXB3Zb0t$wtHs3bzl$sj=Q+GdZ#*G7$U(# zOFJ0G~cb1ouV$G!k$YJkZ z6Ar4G*3pTk_mt}2htC}6?R=5)en6V0(wcnRy9S9ev6!FVG@3X)IE_-gY$ud9;M6`6 z|BoF-dG899(LJSs%!|S47kY)l3(lgLm`Fv^RbStG7xv-X z_=YQfAR@9=WBE|GrrYvwnwWXHE@r)?i;ts(GSUSVr)gelpZOwkNJgAJFY02Wa;d*p z=-Vkfy2P@hG{OT(W7;3CpB_%VH9oHNkSSy;h?^_)!nQC{RH__%nJ(tWZL!+I;#ym$ z=740+3%|I6d|qazKVJla40!@JcGKdRQe3Wp0FKs`I^*NwZ4|MyXu3w*NFnQkNin)G zBud$Vy!QS88KZ1reYYp~<4-nu|D8X!D7)J^A#fK#IEZzjNy zxa3r$-#K5gBuG!(XY+~OqWW@SyJ>mjhAUQmB%Zhb{PUUwvunC)_m!S>ZF$xi*2kg` zzF-`(369lxs^g!}3nD0t(Y_3cI zdHAS<73qA@?)-tNtZC7cg)gICf0^oE$C`3JV*BTJzVC-(R6-uFT@5BzW0O}7#F4#! zedI;nAzOWFYyo4FmmXNM0vt`fP6)RFD2vJ@cU-5&97K`Nv~eLYSO_I z@Y=Kc=rH{jhQ-T8&Oz}IPj>JoXSrfKIF~1hbR2euS^?qD6``Dh zM4cEYg>@!e3u-t&WAr;e74FH9>eY!CNXW>DQH)~PZVg8AtWt1LMDEdes#BYk>mJMn zUyk5-sc1G>%B2#At(SC4o2`N@KdP0Z*_eMXP-0leVAbGzWon(Kv(*@gUE`CUw@{PFzHGx@*n*KZW-l(wB5{#DOYrm5$C@&+ui zSo{54B+XK9MmxN|-G7(W?09cO#OH|GX5>SVYWjy}X1589%CR5QoQvSzhin!iO+GzR z*_cxi$IJ>0#C#*_LU=z_V76HKB_=zi5?t+BJF^@1$ybWFF!Zple=y-jQR?&3UHCW& z7!(8ohyr7FoxXO}5v```d;Qh|H&gC)KGim&ya;3YeGpWd*}M4V)&CufDSxvQx*pTo zpnB%RL$GI9X&lnP9HBNql-6yz?EIx*&z{EZD<>8$J`If$vSQBgJU_!XB2-c0<+*bW z#-tZni{MkyY~@d@BEJ6~oE%IF>z;2@>Wl2too()p^nG!$J^OR3b)PWN(t&)A8ei|d zc8P8|rkLycOmT#O=Ra!!`w9B_wkWxp=L-im{xCuL-PNI6Ry|A;UYm8^H0%W7i}c))+j$Vg<= z44c88iWdxmCA3a;?+dofR&94YvaoNO;iY0u020-?HexLhYfhm#LK zS7x@|PgYSjkL$hfgiN{6v3Snuy`)&0l^R)7&wJ#8knaAaBRyh^XwgZ^{@Gn6lP?_X z{X>z6HQ{enhlPC=BVWS_{}t6YVc>>$hD<0FAoy}G5!FuqaZAfh%DEdeH(2;k;BFX_ z22QJfS!3|-pMowA+Q?RihG)aOu|9G?&_SoIjtyV=0Xfy`JZ|>S6-=SU(dnWOW>tsJ z5zj=Jd`pAsR6l3v|MS4(Fu>r46R^vwz(RALK^L{W?`>Nwpp9-F1 z9Hmu$m-(+bzb z2o|N$mv;Le%aDJQ{PVN*&7EYg$9G@S2e@wh+8#17SyI|P<_-0jeSwTQIq^n7sQqt3 zPXtnBjRFJgD-|n{C%Irqu2N9}~*v9LW9fZI??o;h-ZDFAQ^ zo4Wfk6)6shGm`f-gi_V0m!@o<38v}VN?#CNd{Y%j^ zO|4m29^Mp3Bz1HRCr8%}eGe__xh+||w`81G25`uvPHq}0Abz}+U_?CIpLEOdIo|Vj z868}DjxDmJIKJGUqYhpr_!ljo;g5zN35nopMO{`iZH9`_2QU(oU1Rk?QBX6=Wly2V zUf}N=)vL~;21>Sp1r4vUaCGgy(@Fii{#qJCIkRGJYZJpO9fIY`4Z`8BoGXXJ_Z&r- zo{6So7)AD5YkJz4dRU0}IS~aS0PfZqlI8!rq@J^UB9&Y= zT-UgAEKpLIsA90UE2OYI13a*&+~(P|hPf=I#FTY#28+coM zR5Y`>aK*38Er8o7P`?_ z43gd=qRz{WI(P4F8MssWbUdHSvg}SJD9>XAU&RlYriGXYd-P_3G>OwN+alvPo@x^r ze@|<@n}H6P0*s!t1sE}ivH2hMgF~w(rmG`4KUZMJgu9d;;%{Z0yhd$}qaG$vimCOr zKveOUp(nz?y(R4C zu;8Muc4hPar>YTcp67wgOWPXt-X8iQe$jlmh(X)Y*L*upXA^f4uIrGFS2-*dBM?Ni zJdFw09tPS6eBJs+RK4;ptuOxp>-uCZopvz|7=}V?1kLFlSU!XOmz68%_Z*z=Pdj!C zj%2#54Zo0FT49rtV*R_=jg>vYh&}XzZud^#pDvd6MJxtl+y&EUBno_6r@l&@IE3Mj zFo<2Z0imv>kb6@dVF_ji33^{P$O$FOb}tDirCmPsIMH%*afnJaEDdcn-7Yqh&<-dF zQuriYdOJ4zY}okTm4{jW5$%CrrTyG8M_SJm>m2H1T9eCUPsNVY-l*p9qy}xMM6^j7 z<((Q173+{SoHF>v^aO|tb8D9P14x@`38(TFL}YsWT~zL4v7VRM8#O!FY9=d(xs)wp zB;I&T8-Mjd8(8}-uz(+_f|1jldF4}?Ce0yPX6@15F5#D=6dR`4LLJ_V?AEmArb?x7 zwY#D`>I}V4w%cFdG7G7|Av$^}V`#_%esBd1LDRb}AMzLBGzW6ltFw^6+&z3h7iNT7 zNIb-|JlTljPr27>7(Tq>+*<$O6+c5%r&%RS2+Xc*CMLJ^*C*Q?RQ#2LwwjPq?g5<) zzpyVlh6D$+Cd9uN(;Rr8KTPAdApC4U8L_)IGGA2uVW*f#YvhALx$mjx*Tcm$_t_V8 z2dfDm{my-Vq>IrZZrV?X)QfTt(ix5(cKudak6Yz6{2E3?$2;@w1rgWd&pbxy8U+w` zFM}z|8q4)Y*!-pw+{{CczG2vF<#~3!XU$nOFCVL(O(~~ee)p3u53RL8Do5y@LL;M( zE7?njg~Qwb>Hc@`ZL=qY!{b=S(JLj(CdZ~ zvr8bKNyMbR4(}O$(U0c}K6uFDmoJIY97XH31%R&yOBsG(Yp766_Bis<*g{=_9hJZg zENr7b%1A?shuCrxrks7}0XQnbTOp*(5o@*1-)>W#b?G!6Y{{2Q?_PaL_CUTaB+Bu3 zjNI$(10U9M1py@`CH~D5-?o3pKQR3>QByq^as8ku?4#`nsnb$VE2G2L0GjLC zTSM388w0R9F^l@&GHO!0elT02+Ew5c?L}Vw9HKVgrqD}-MD~MOC{YZ`i4G8?Ey+0HRg2)9>V2n`fjs43Zi;b8h{3Xt)nll>`|MBR#k25uh!SQxcieUKjD3z2&ynr|| z(#k;%g{s$_m^DpV$$Z4Zrq0We8~F9BxVe>8{)sqOoF@$H$<@AOBXttwa9D&eDlTC} zYjBN7CtD!Z%%nSX(^R)Ndi1je7Q%QuiNv>RUM9(co*38q)34X{y*ynkHU?TcdffzX@^$u5)2p|Y6Qtv4;OK9Z3Hs9S=1xkB zv&8Ipi5Pw5TzFo?c3d$^vwQT7{cw0WY-yU}!G@U2-($T}lgog&?$vbMYPA~<*sQg+9o zBS-5;HN1OM-X0|${KiI}cHM-PW5nsIc<~*%*vvfF-vBmG5&>*?lm6t&UIn7nT7cWA zWSuB}GBoeZGUoBb1mm8!U6Xd5LpnQPxY|k5Iz!eSv+eczL+aHADZR4!_cArtx1!0t z`FGuaKBd27!|}^)HAIlLKGKtf@h&I%bBA->!-kVay5m<`EA1YEY{oN<0S61oMrAfb z-1VZ@&?#+>Vw5bbPhthPR#B#2rKRujHm|`KZ28v`#RdC+A;)a(OuUYEmh}vX%pn^U z@xyxG()<6!6(I!#Zy@wkl+ z@1x_qA{s3Co$k^ZY31UQFa(PWFMp{0^DW(Ns&K*(o#{Wwh=whC9o8<3!ddHw zsEv3;0>v0w}YJ2Lt`xs*T6`n4A$OzNXc6C)%()#l=ufZZMqOJEj ziRtO@CrFVKy?76Q?H_~FR4+A@aQ-tfR;nM>@SuRBhPAdAVqD^SLNK;0a1Oc^xUyCn z2kFsgEK!Q{`}&rZmsIl~KGDSS;t6V+*>_+4H$u#ce_veVlLo5##QbE)Yg@bpHSQ%-swfnJngpVHxIEO5t zF94SFO5R~2h5SmX9;Y-Yy(;LE$TK#fr@wMb;t0>y~n6H z?DyJ?uEW(zZ;9Ba9H#4iL|{03o72FYXG10^&y#}1PCFsJ+iJZu+;NMZxJ;MlS&hE( zW~8Wf;OZ)ph-Xkj(9Wp8O07^eKd60P{ev=Ywy*c@5Xyy(dO1>*sgMPqSu{3u``!30 z$2M;!c={N3sukA0cvV{_=U3vvoGdYJ&@_#0lW;O{g~mV_yq5BlOyV_&m~9R{p5&wW zk#al#d#jR8sbLh_-T0lKDeJkq3jL|uH*7e~E`$?OyIJ3iMJDdQ*>FC!DQQzeC#+UR z5wSW{ob$NKzCbMtvq<}W#HHJoTeZ|3#y6>Ls)+*q&(HQ7Mn67z7#=h?ZeL51;k`rf z?p{jvKgHXCC$DRlG=g6`>C6ehk%jyU$D8D5TNLB`TdzN}Zc*Gx%hKI8pT1cO%(c=P z$%CIG+QonJ>Cu9BJMW{ywL!XXy75YK|A4tmNDJR~9Y^7dwp)dyT43b=nu}!~hQu?u{ z>eQ#6y<-a3ht1>za|XwdDd3t1D34{$NjMKes>gh;xsIuYNgI7AQBAV2_|C-ote`>7 z1!!wNeAz}sh>1~Wg#*VkHcKN>6Ce$^!Yaq?AnojIBa4i;(9@%O;n!W5JL!^4g&eeX z_js)cv-W=M-Ymqyj%{5Q{~eC!aETsOcbaf}vXSqK4(W7xk5LETY9-xwI67!Yw6F*M z@pM%#kJSI5*BtdDO_&UE)AGk_j-~Y6R5dKNgR>=YB!|ZRgk zpe9Ep{;WWL5cdO(47#iZXjea(G>@Z6j|5p@2BPTFy5Av{fFpr(0y!fLoE#g2H79xZETjy>E1K)UC2ny1t|aj#IwJPa9gu zXL~jf*{h_ipZ*G;4`9pe*F*;eXl zCZR?fso-&MuX=Ul(z0MkgdoRgQV}m%%(3uR;YisnQs(cGYkk#fVcGaWO_|suP4bjH z7gd#n`CCa%rAPE3tGglNSXAL~AXrI5v82vVDmAVgY5Xn@jj;LjE0(brQRcew=;Yte z2c14&4ISj=@R~h(o)vfn6<|DiUVbE%hozi;KkgC~Fyvo#TLu9W$i#NPIMvggd0e)M zk~4PeC6ge=_-DJY$E9jq=6c$moRPM{9YyR>O=T=>v^UGm6p}1BU&mN%*+Tq3Bc%;% z{Y9grkD*Sz4I%5g1!UT7Xo3{XNN-H?^xs1sqxy->`X>Fa?iBOwF`?SUTAkdQ1+h1Y zV*ZSlV{Lx@ONz+Va!}*J`EpZHRxW@XWeGJ2N`*kaoWu@u?h1UfpcYraZuDltjmO#@JOl_jUVVI9j8R?FIvFyFf3efH$Xo$H?-wgW9+x><#tWzi zIMIp4rK5t-tnDZL^AzHh(qX%bHToZ_jSP$`%n8uZ`{kwEA*F<4t&wz&OX+TNkfA6- zN}aF+i3*;|$n=4s#9rzcMhF(4SR3 z37)(ZtT$UJRpuzTeci3Aoig3Qt$d?1@(zFqc?oJg;1=iIt_{Gk>@PRFdvdf^3TezI zqx6xr8|wZv_e+GA!R0m$=@@nEP*J)8@k5RODpvvU!u7G_fh2ND3E*7<`D!yg>o0Yg zelNqp9qPbFq9T+P?+jlyuNqG@acl}?tQ^y~_GKc<4xGG7O~b9(IM*7f4s3(-+Azuf zgwuoanEjRf4AZ_$8G0c*`XbFj%4Y}TzdPa&vo!#1ml-!<`g;aNe6$;@2p%11>WIJU zcyizfndd;~R)s|Z0m0AzDfIMxd~UE_=Rs=39Q*t6dB_7>6(z^f?cvM`8uz?%JRk9avjTZ;Y3o>QlxY#7IyL>w zHdUs3JPih8n(+v{A5b}z&P&+zkI6s#1hV<53OASoE9cmyqy$_yo>dC5?XhWRT4PX{?-O}Rx{Cq+A zUWO#|MYOE5jA9gSLsyqJ8H52*DayzMJU{;a>avK(4w9dr|70&R8Pa)aa_o;!$7vC2 z5DWipjO)~-6T1GH2p@j`7+q8eYvX9@Z{a_D z81SRfO-rFPEW>}s3nf=owr(KpD#sGo3Udhw3Cmz(N{kv7Gz%Dz~WZ#{jxnYAVgJ3hgHG4 zc%Pr{=h5PCYuqA#RB1_a^lv2yNSed!Pg#N!P3zzR^xumt!89woYAcon5kCgQbZ}Etxc;K%|rl&;_IbfDaXdf z&=s##fhU0gRxV$^UaUs0`9c2#+ttUDRqEK5ff5wKs2SD$lQ+`4 zOXKs~jllC~xC}G$M3@dro>9e&0cuFN zLb>bikWMpT(8Bz;bH&Qca6YKPJBH=SHYXrsG)1PgOAR-YDK=n)qqS~)kYTVkoX)kb z?H*YSesSOC=uEn?R`q}lTzTb=7>?|#@VjjP01VE9okVq<11W_$?Oh8@7Rlwi(ZB!S z*B77K`!g8{NO@+W#>ra?1R7lY+xWqzpip>zuw{Z~0@g-KAAC2o2Up3(Wy^G-@*z<1 z^_4`Z4U5DL!W1rqJtkSwVY!%4pyx=6kV{mWQ*NE#0MZ-cuhhX@E?5&|eHe8RNmQE9 zQwF0%YVcLo*7~VlPyo2;w0MtG&^J2^&dr*>O%S_pezpRD{8g$dO_hA7-C@dVjiZ8zQvd!_RXxeL4II~VtA1y5QI5#)`O-`i4ULby6^T6YMzzO7b_I?xwckCw564{g%n#J*T2r~s=B zhZG`F<-FIpgRFN0@3-jVE93xgN*Jj z11jquPkFmDq-Sg`p3g?jf{%2?FVqpBd}5CJPIm6N0gNQrppf$hp{2XO|R^Z9M*;-R8HE-9V|Z@?))s)vcx zz7oOFxLaQ~%}g3fxr}WJ)vYC+&32Zbv(FlhZ-RRd%Pzpyf7xWbH9oitHD7`^qib`U z4?v3H;~~z%B;hHLdqYY}$|Q(QRu-Z{Qco)iS}g*InRgjzwQU$Sbt)Vt!$UHWt3IZe?+2Zdew!XQyurIoH9 zKjeX2iGN5sYUIK&^*@_7l{6R5b^GR8_ zmMbFG;u84c6f@*SNgnVcCdY6ZltK_j2udD&b*XyK=-3O8s_n>M zP#SHhQCnpXc=VoDHUIXqElQxn`zwX45{Jn&>-sU0roT^wk`R;+Pr|!b)X5@$uI*mv zygD3I$zm)FK%#=aEf)6lR~FFtAaLsuVcpy6`@}fkFq!m{f=IKRM3`V~q@m5FcL`hn z;g>I2`WdNXTmq%lpyD(}FhIY7*Q)t{yHX)_0~rkf3bb~|>02Q4fv2u6Sqvbxt_ALd z{wPjumtA}&QBg?8E!3If4_o@jv}Xn3$yRJVhLlsLm6f=2SlB6V?a~}AJtmzmCJ6r`Gu@rV0@yV?*3-8L5Gzb}#aTI!HL`WC z>LnIq+kFXHq!uKu9BqJb=PsmI;iS8s9z1Mt^BAc!5eRH(t3`v-uv6E?n^#W0g8m1& zC|R9^`Aio7+O^N0-+10~yim5BtPa}Fa}xc}DMx!glo@Ui`@(Qi4%)-Sal$L0peffR z3dAP>X7<82g6I*Fhzrljj@*DnrANl|Lx2pcf&GvuX9PRa1FRaae?xR_% z8_yRVF)^VJbVq&l;Qb8AbiRwgfIt<`UUJ6{K}L|{yh5SiCEjj9fCzix18P?7KHzM#wxS%MP~Q87uyH~yx)2{ zkL2^yUjd#ZrSfp6xxTaJ(Jzi-h|FGQ2_gyAy7e!2pnF`O|q<@$ETxyUA+iU523-K0T&C9`7+(N=hdj)mtk`9Bq?i5C_ZO@XBsaatrcN)$j7jbFFh7(pHfqgRrK)2^1V z+E)#C%Qkd17_|L8ZH9%|tyr52r)!&sbrtz|US>3Tqu%hXq}#ui)cpqOIz^X6<0s?q zn0=3NcgB7tf;NrfDd>wl>)!e)K^)Kj84js7i$y2(>F0{KEW?byokalyf=c*=@>ZGQ z6{uIlw9oI;&szWbcnt(e1$O9K-mA-ZPd-C-#XSyr$jN8O&W@j6ReiHHuKHGdr3%Pr zPK(U|qdW%ThSr*z3D!%BuFiG@g#ZhO00Z#B&<7qdXSExs96CiElt~U;>C^)otFAJ;5rT(v z1>d{7ZwT4yi=*voAQx4=@p(LC(4vwBF+J*m@M_f4I|BuQ4?^|0yW#6f&^aKF?#YbB}xD4rFAlsH!a0`Tr zivi_oqAB;V*y{iLQ%Dp*Ahb00BD^9~q~});$Ac5zNw;&c+c?|B%3mnlk(a}BX zm`Q0zFJAnCtxeotJ)Ibc-dwLhlssnH;{GPW%$(6>V}d zUfbu=iJRx!Q{>Zm)+Kxz{o$d|fdmUXrwM-wp#eK*aI87h-x7&LfR=?M=Rf`$O2~BH z;5ZFS_pf8^;=kg8Wr`&iEP63_8vbj#)VCWBEp z^|Oth5pcBuUBJ~Z?2G*PobN7&+2D`~&t{hZpllBkBOAPH4yoUQ`9OB~~K~Y~0bX8*bCWPP3e0axt8htNo zu^sCHL&)IYlwSu<6t(Yptz=?(x~@s^4FDgr1Z%1e}ax1z+-x)(?2jY*Ua0+EmdBvj6{ zL+(e84MMMvZ|zF@t=Uj67kjgeg=y+6 zCpon0mO}TYhr`cJ`iUCJ3R^3V^Y40j!sx>T`%7RDv23AC?R`1KpoVY||LV3`FGw`H;;yEoKWZGn;sWaDy@yl5SXtSCLwwK&$u zRK^b%A`HeSXpZIGCBnG_C#CtI{!v2(Z@10iMwMglt%dseBcY~%oMf!bvCEHOD0evN zJz;F_d;}U#yD&n?w{=EQC*|h^wu8lt`>QqbC%!An&{N}-H+oZ+k}=TXd)m>}>q{|O z_0i|)$@KsQ!6ZV z(xm6+YzE-5G={SJhJl`8M@z;BxUM`&&6f>dn6GDZ9SM@_Vg;R=?|QR2fFVrTuHKJ$ z(_#J%s^#ivfW~6E=apH3R-fwW=_z}N5(%Ul?}c?;y1xV=1nt30#=zGrxNh7-d%Ol| zF2jsZ?kbDYO5MPG)k%1>we&pL!jV+W>F3)JA_^YGVW#v?mLGz-)Jb^l{f>-R%>L?| z|8oUrrfOx^(^ohS0UG}in@^sT_Ej2|Frp0ms&U)Ov9CgD-C(zM@XMAz<*Rn_(Aw^4 z>53mCBO2l+@fl4W52)~MYdL$}!T_vti6-`HE-^*nle z1CvvQjQ~~p$RstZ_(OOHbhhoizt*C0>u>Ws{jfrB*{N*IpF||5R z7teE!mR^iN2bMmJp)by9d*Fhq+)U=*&ifuWRoOM~++z0Wk}*fS9=i!A;;~QHO_Y+4 z=m+wAHOz0_cexUe>e%Th>%YbGlVd#eh86Vf82`1PfcpJ4vG;5YK+h(WXF29ceB{xn zsi|>aOtCr{_CokWP2zYm8gExW%lwSISMO7)CZpT-J5u}6z{Coh=vXAm*Jbs%zpr{3 zG6?1m4QHNTk5-GlwI)bX_OaS_hq7%eAECo5I^rJsz=?Y;kXvV z_1Q#q{e6UV4XbxlgY(hLya?JGkI?$G*J5tGMhdcnfOvRV|FSM{1RWhKN(e;zd>)-( z7ZQkmZ&2m^b5mDxQ3~PTJ~j4&Zqr`R8sDG-m6+xDQj`GGFYC2@v}dpqdwUFPJ3!dP z8^XfL)1176b+R3!j34{mIHN6AP;`MdLk@dVn&*6nqMy~tsSRvLKD6#_j63#9>~Cr1 z!30v3XDU!MTaEjADU~*?izVC>#$%f*EYE9hxiK;N#N4tPlKjEauTf)qUqGR9ofnUL)No_7M*@H{6g1`kfogSau>fz7~A^1zq0eR($>oK``1l?bX5iKuW9lMoueS~Pag2dBBZ%~!XXua@>mfD6Sk5>#?X7mn zu|Kdzw(UsI%_%I*j zbNz)6um0PdW^E27TiQ%58@NG!G#ko#)91*UN!&fQlu$VisacfZsTZk?E>Ref>}nk*Y0a#M^!I4Y<7=23oLBMYaP`u)*^Pzfxz8{@ARDFE9u z9{JE!8Wh~K;fj_=>wa~=6+*(;rW21}=D*xp`#6aUi%IH@Jz#N4`}<Dn#zp9NQ3P@8q3K+FtlYdSi+I%i2IiQKy^(3U!+aHYNGghlXj7}=m7#_e z0x^~2>Qy}SOr3fqLZN$L{@=<^5yv@fP!k4&lr6`q_pR<1=W7*}mD@XAe1FO0DE*sf z8-e_CY$g#7{+f6#J~lr)7wQrYvz^yV?-j7&*Ru<6;qck4uioYZZ!mrHOP5aVnD08b zD@5kLvb_b=RDYT62(LKk(?y^GOlkbx{j`HnHg=G|SodXO;K-;ud-NFGN5_gEdY@uS zNxJ~cz4H6bv=j3t@|Qx0wZ_#6l<Tw%nT!Cl_{rjpZ zEUti=b+fK^}~sekdOsI2DvYcH1|GBV%Rdj9c7ry*27fvD4`TdBn7I# z?=s{FkXKv+QN$f++VD@a|8gzsf%D3QUBd>of~cz*$QLZW_x3*c zJ|C~&CTVqI^+^D>RLABDy{EkXi0G|8BAy&-zpyW1ZK28}MM!Zza5ThSaIna&P3B0I#t; zUd4x&A*Tx3?hHL)A8ubNZ8(cSca1`y?s+JSgy+_fdVECt)k7V1x%t60vYj{&3E7a4 zO+{9qm}$F!$OJ;sh(>FhFc_^N#ba@2D~hivtyyHS{`83F(Z{gBER71UeR~0iDW++^ za|uybp2Ig$n(=9Ak>I5=_@Yx)l6hz0uoi$u&xl44>Q#|Phc_s`o*k{d7IEZ=W|ivE z$s-Y!MDcTRlVc)4vU#(>t#%T{jLZ=Cpak82Bba1Jz0dK^5KrMZ&o zWT7l8d9-|4{QTtaXqgF?krqcJ9bd?hjY15I#1`xle1;T!0D7J?fPFHJ1^jj2u+Ci@ zGxAe^>wEMB-7;eGzW0?)-06ARdZ`o(oL!ZI8mU*1Tl+>j|Ha%T%#(HK>q!^uhYnG@ z>3Zg;zdp_uqk^rY)v)IzZuMTYfCG`=E`T4AG&hyOsgl%U{^=O>o%5ix5piB7heo|# zoxArP{a6{}C(2CZ(B)&0kZ61U(=Y%h%zyj#HGJrli1UusOakYDrZd{7P-iZs&}h>K zUMe%Rv}(-GwL_K;1v^NI$<0`gEYRzxTlSdY+&yPa|wB*-k3Ow9b9~TWUxC*<9A+` z#|Z?GBVGy{uzXl>w@BLYdSaq;8fmkYUunGAJ5}ceU3)hM-#;F4@H_ZYJEYr|{7m0% z(xo$b8uFSRox4$!_GI+r<%%aYPA_h^yw|G?;ZwfD`Sf1M@%fc)9W^!oJ5GyS<8hoS z2L>E*utaSQbwrgqiK^Tu&!P={M@djj<)V`DlaO>CW+%KO_0`A{#y5GJqqH2yOm{jt+>W6sXv+8j%6kVUDPR4uGwAD(-IP&JV>mxUuMLPze2>dc z(OOicvbO=xO6R~y;X_~aT^Sh+SWrAvobp$i-TiG03Zqij1AyVL3@7xi+Gym&so<*oo9|ES&PgcTCmZ*Ou`ug+-mqvCmH`vSyY z3QtCJ)<_``_BI|kNx#qhJxQ7wLw**b?*D7ZyZeJ7cjMQm=0qLC6=A)$DhLGsBAA+qnd`Fq7W9lk<^(kfbwhE@PysN!kpPL6QDn)g0> z_y*`|ldYL16B38Y@r^BaasO2goGx5VP0LK_Xch_3wP;JnA=f8rsK83vK)4u$)*C~V zzffrJZA98ne%W@5iZsvo!5m}G9mS#-dXnt07-_Z&2%*5JK@xxytA-qWS%tj7_Ls~| z1X>~vg;nUuU2yZZ`Si$nb&$a4=wB0w@I>RZAn)?Qy{CtFduGesTulyNw{$OGTCs{e zgeIR0fh4qsfpc$%eh{@!Hd<(D7pSu2H>c|QOO2$?>`}AOhSjE}WAFk59Z{tbX?`~; z3{P6Mb_EE%0cQanoi*mz+ef5Ae&EfP;Xl* zj#E87+69s@;pw+vf)`Uk*`KqgGN3~M2i=%oj4>?<8lX6&!W-&2@MoZW1tX>d9eWrK z5iM8vmF!PA9B4JMRZFRumOTGmrP2{-r-ibM3_13Kh;fyie-(THT6~2T)laTX&ps1MyRUbZZK4s~p5^|xTR3qru4>v>EQI*d zw@unxdP!MBHi_fF98GSn9)X?B!50k_c^URN$8(z{V%_ol$AxC=p#hYs&U1IJSniu= z=r=vSwIDLdMqu7=f-3Z0s=&3GliPq6>hL041EJ@O0wgEDi_Na6_mGv`0tzJ2Y}m4r z=GY^$?J!FuumAGDCOlR)z>a9CIy6!q+W4n60O@!t*`0>f6)TX?QUmb1xIWWwG}H|n!9 z#ro+J5X1rc)hw`jkq)%_RUoe~N`8c1`17EIL$}s0xR%U$4`cm?ZWGmCVCVqWA#p1( zI^X#{2t^v*_rsnup-oh2o6vLNrhNOe^=LV!d^%9CSHpO2ZAb20kWay6VCU5B=q4|RKw zAgqN@W6QXde}<+}=w($PwR%8}_DwZTOW#-=Wg?OV0{S2%szHB1$>EkkEq9KO+ez5KlypZ!cg6D z6FSp;w$S)0vAj{%XWXEGE9Gnoe1dut$1LMpw;;I9wjY& z6!Cj^`z#yNPb?M^1hhYQ-)h2xZ{#>!js5ZF4pnmqiCly4Ddm}|u0XvJ${0n3Kgqwv z3++KKp-lw`A>z!&8cJp;dGkt7*yLPFfV?L#!L6|$#LL)CXJak&Rvd3KMBUmDJq5eH z>Ab9FZfQ9;buNkD_0MX%v7o;IBu8l92T9S6+qM&Q>(56CPt8OUU0X%ce=jbrkGOI_ zmF@OzD#ZUq)p^Hb+5Z3kA|xXtE6FB8C<$eg)j)R1N=XA1N%kHYiKJ2?rN}B#WMmbU z6h$Q~DU`jkey_v*`Tp^HJlx%Pb#qC+yvOVHd=sw)=}<467Bl_!S~;zMU?6i} znhH@e+Hr1$=gc(I7mmNufm;^~kX}AHX?FYg+A zpS@4)($H6$$9hjL&2tX9uS)NIm%m}h>00c21-^4^T%(rCfPnwh;7w4%w&YmKozp*lE!#@58=claAVS z%bS{f{^#d@zqHgn0*7ULkuxJM_AARQX|LoL***8AB{d#0&a}(qc`Vhg+i(B9{6gcQ zxHjXuTXdbjdyZU>imHCRu^%x;wXkR5EF|J{yJXgk3vt3ri7Kpy zVnxzqzsDD3T%AAjy$jb>N<{nD*tB!&OAo^q5x@Uh;A`26sNG+3W?j{-JvUTsZKf)7 zVuRtP+Hw`e=IO~Ok_l^vGgTvQdu18n{N|c`JDjoPl+BXH}!0j#S^Xr{qg6& z6YFr>Zt;ok5EjSfrN#Eg$5_4xEc;cc9I1P;^x@uG-_{8^t2Dyq2tp^SDn`DidsU8o zdxE_QN8fX$DH5!KniA?e%r#R2a-5R=)vZg^mIaKq)tC5`*@^a>U=~nE8EX?duw1$B za#O!)MISgbXQT`Iw5P>BIjBprK3`Zbwf16V`+$^5Jk_*H$Q+7D?mu3!eZ%7N?_U_$ z{H^&jLcKBE&>sBzn!0+$M;R@)d<}`4GfWg8Ud2AO39kyJ2gZU)f;8s$I97x;`QLsv zj6U{r?rJ+2zE{qxWE%^)*0375WO!z9>9W~8Qv8+DgG_hfL*L6C z-x50`PT`Mo9cU30@Az@bunn^G@{B?rj}7jV4J)2}1J?V^vReE-q$iPBYuJxfW}J z$L5a=&U?D!7RjhHU?7qz$^mUhqhx#uYDdnZUWpuP-xP&5f3AC(d zG$K9?B7KyY8Owc&2{dozEufw+OTdgJ8XGtl)rEiKHUx8(ifvL1Dw zEXONkVDB?@85f>RY$DzU=f+uqMx6E8U8mSD&ap|~4IO>eTd3U(1v> zW|O5ub@}m}n_=rwc-0y&8J%%v9yh<#y=R>F3<}lp9rJe}e?$7WwA=)J?`qG%R>`FL zG}>E0ABuNrwO!`B4Fi7K8DIagI&NRjSM5VAVuFdz)@ypZ4{A%^3^!G>Vtp6v=5=uK z0=Hl4+Kw540NDVub?cz$x)o6k&^z-ztd)~@yHi<9f%n~5i-9b}8zFkg~9 z00}N%+4+%r6lX`_G|xn>48rVmyN4$!*#L^*kFQ$1(5ze-XI#5crIx#)4Id(E=iky* zuznI}=iFccjmjKD#T+%%^xJTxhca(q$loabFWk3)e{Htfru2W1vlC01c4iAX?&DRb zc`A*JhqteGwJQ4D1sf`gM6X5`>`T7lxk{4iLUl|w*|$`*>-~NAh}Vq|85MG=8n@&X z%*gH7!Hr0!4wk4#-|~^8mLm2}{yMsScD+~F)lIYP&+$2xyO*y>4oVF+2n6WF@@{8% z(bYq;9u$SS&)A!^jy?E1BoD0_6ql7gvrEnB`=2h^+1oeqald+W9+Ul`Nt@72o2_A!33;B(W`EdpZzE>D!eS}J`yxC?SZGThbe zZu?EI=l7p=5R;@@t_Yv7y38=im$JoH3pM~!vCc#+?dxepty*CN4QAIVuhR&>2n96w z;m6t+l)(Eky8BZjld+brCm`4@zoK+uTA)H{^O8ubM6rg%_cuqktzZoGUY!Vo;2ut- z(dphS@5Rjgh3WpcKCTM87J(W%{pLX^|JR!ii5zK+lYRy$x~qD6(kBnZKRH z#e`~C1a`ctp{Wd9DkNTX_-VufgjD554IT?U8~LX$#|!GH|2E2+)-lmM-#a_$4D1Dt$Q7( zC5@_@^@!!&Af4R!W~>Z8m|Zfv#w1so;U^UjzT$3;HvQ1|SOPse9tD?#qI(;ua7=`w zYdc!U0JhWUWBxc;Y&PsyUsNEe@!!s}%tO!4==%SSSvOn}EJzSjRtT=1drI9Dm3m)S z;9mQ0BKeJ#pijwh%RQ^o-34TwSEV1hsvneKkrP(8qd}te3hWix_%`q9XVzAK8FrO3 zSKLY$sO+Sm61CVE>zAi}!$QV5b2qc}p;A^fFlKTFI5{}Fk#5yvSv@+}!+qbeE}A&J z>o2Xp!u7kznv#@~5`mJu!@}$s=@gx3bRk(b05e+r!8{LPg=Hae1_=R*oSWztg&D(? zs}8Yyv7d&r@|N$of$x-g_nn)zoR<{j@ZrSxFAmE;acQAOHMu51NE^J4V}hYl$1nNt zuo0LmN)xCxYYnXSYBNc^8);W!le^*ix=s5gZ)D$0OWP=?OmR}Y!D!o44B`c6J+4?p zr^vabg$rWzMI9$ec2H^bij-K%_KXs6e-1T+m+zO8iYsASM1PQCGVZo^jVxcy7>Ii-Eu)^nU!YooKuiC)~V zYf_7zWysp~OkU|iX2hf_zW>%Gyi0pH(Y$ASBIayg3eWv)eo%OL4ax;-e*f$dfR*-N zJ-9IKU{$s$@%slr$5- zG^>-<;e7dHFCHof*wcUE*e6g70=i^w3Vi5bhPP|>$c7WX#JUstFRqSjzBiJ6L8Hk^ zoa;Y+{2){)sAAVD(?@5a{Q>L?w*RIVkAhvd7?$3+-hIOUdFSbCd&CvcM0zE-@fm9R>r@5V>2}x>HyyFA$FfAY7@^f*tn_uPuc0N^`mt7j zwc_^bl#>9TqC$W(#2%am|M_gtWi!r>L#uI*$yWLNUp^E%YDd_wYAmzx8@#XSxE^!M_9Lg7J|>4BCEcap=Fi zYj9;0^~AoooqxLpgtH1LE?DtAAZ0)W#4b*8%uv02z<9e{o1N?48;Zoli-v+?$Qa`nky9>%+Vp`u)6Vw?K25HP=w{>rDgcmxJr%ECLbo!~bn4gpbhNBl zUVcC9-9uYBz{uUDbl~Ismi8*{1v#W~%u1>Au;|xdpCwtnf*ZdLS5t=0nc;`=`={o= z;xsXPU@B&RhR4g7?u2n@FM`?0p~dLQS?R&$>MOpftE@+;Yk1QZn(e`709Ncw@vY9E z@y*Dlq=xtU-J}`RWI6)ZIN6ntzPlL(8~Zb3S?i)T)9O&uAaV4cf4SrQ=0vA4LLY?3 zx23}GBA@aZ8nm2iaC4KlFJz==)6EutKD^;cgE>bXd)o4jlkMG~_pO<-iZ0fFXCh|O zI^TTp%Lk_Y2d`_29C5ttS7|3F(sV5O(wW@P1H1h;{+cWnqECHK5009uF8^GvDZd-^ z^9g%yDzUxWF3`7OfVXmdW0#xO+&Jgx-{rS8Wy?GE|7L4{I(I_)!UX3b5GtUR$BK&S zTGUETPNq%Ux3arqy9dOKt{zb?jI8&4S>7Y^Mcs6f`VRUawNs%y4_zJ!(b!R%cgVGo zgk7%F=gu5s&(V$GA1B+@=acGGTg1zwkHa`9)jNhG|J*;kEvzl7f};SmCFaMnnTPK5 zFpj<4e@`avZA#wG_W9~yC8+da3E>1E1hTF(fuZut^Igj;F#&q3R0H%n!ryl9_GmbU zlyp^|Hm8I>%fPZXFh;(}S`34);lbYR>&rE@Z0uCc%-Bk6cU6(A)u3Wd8!b;r=2X>s z7HnHAzKuO|X%WVe!#>lwEU0$=RQ^NsWb2uGnF{E7xWnu z7Oj+J-lS0#$Lj$pu8h>n#Y*XP8kg$1-&b)DOBV@LeT{3%RKyNFYpOz>a*kum@2(43 zu0_6F{BzO@TH<(;{ZtLc(lT|=hGtxP%f>5bPk~?JJ37FDPPv5FvQ+Y&p=o%Bhj#Ol z|BT7}^w88bl_mZsZlBkpu8xxOcGS$?KQVALHk}w~ZomVK)T%yiV+;;(-9Mr{Im%k7 zHNQm3%&PZRh_Nq(Y@}=jR%C8c#~&}kLy@YHKFck@N89JFG&=maIRkVZN+2q!PPLJ2 zCaS0BltV$r%4{JD)*k&0P=0n~O#7nyt%W`@DQ!-U2rS&>ihFdhBm=aZS4~aYz`s$MC(?M<92j^0 z!8jf8BPt^@%P3`Y!!amsf4Di^rw7-UiUc^0w9$WnI4Cynp-7 z^M7m6J0%_N9@P2ko*e(k;CwRnJ}He5{!@FQ6IIu-*C!R{Hy#*PIvb&se+tOd6~ZW< zalzdK%le%V2$A##$IfW>P5lSz&FQ@Z^kY_46InyqTbuP%H_=mQgw_@1ua0cH@Z)hp z|6_m8^vHTXA%@TlX`~@2CqRhUJo|hlBiaiO?)+#)%^GSi=Dg* zE|gkHR>!lqZuw1HS4O=MgJEgQ6W2Bd#K}%kA1iY3fccS-ji@xMQuc@@;Pw7x`K&#P zrU`-Xk;D?QtYkj{pC$LLJ7i%BKdDJ^1&iU_ENPLDDyD{~FT_!78!Efjj%AfxHhFWt zJ7&Ln^y3B=C)H9ZV(wWk_ zO56l&p11^|Z!f{^?dY)-Z-x?aOWAohl(w(rN0;?d3(c}ixG1FhpAB&gqfQqs0Q4Nt zA>oIqYQ3>;zxb5gDdaQ^3aRYt4qD7Ul(l(Di1_h0@FU2k87tSWh+VP_>gfg`2 zmm_Cml2j-@lP!8gb5rTbnIXv)^d8%d3C<{AB--TswA;9@%Xfm=XXbA#ioSaN?LNQK zBSztWmxwh*B@@rQ1IkMad1CCY+E=F5u6P~pPg0`@0vh8T%dbYkw06UT0fN#ezCpk9 zkSV0UnWbeK0sl&1f&e9thO4FM=v!W_yVXeWWlVl5e~10a{ZH{{NA$&bz&*9wcohPN z48+r==YA({-PDL24dI^&GXq;a(pmS1YFuxRK0fnJXMWaf<|@~D-+k8=)~upqGTx$g z^w=>ds$RNOtY)+ zjoov1(1;cOwY!(N?AV&lYt$URi&IQQM0-uJ<0<6YB%p_~qBTTtHm(ysE6VT#p19ufMq;wEEA176nzomU^39!)i`MGN3Nb9!H2g|ZPw0!LY+bONmH<%yTRpP-16eK8_fH5rJqZ;dL{Ao$z@9Foj_$;~p~ZVE+`tknGX>P+os5W$eb*p5)Mk~j;7yNx0n9TyfiWAub+?|lYd zUT!Be3F8Lk?tU1*y}S@1E9_u{%X75(RxE_{%&-^MgCODb$(>%DOS=!6#!%7*KGuBr zFulxobWkHeOR&eLHA@zT>Al;CY8`;GRPr3;Kw8yCe1A-Bd5wGV*Ul;J7r4reFNdpN1tEuZ!=*eLZ#6I+G)@EsEz?Cc2rDo(s4r~vmZo-i2N;VitKl? zn>l1TBkJKF*T3`vN@i~qS;uhydM@knSL0K!%A?_DDWnfY7laA|Q*^s0%snQ^XF?Q5 z{1I}Fwe)1$hxApC6e^#dw7#k2L49$0xSB9J`kT5!q{~`x%t7s2JK{mil!ssVIJ`;_ z^6zg+jAY(>^~D}m8Zm8Kxgg%_cp-$rwANP(p-Z`TXpOVae$7LY#Au<)jQuq08+xbm#P7im4pdQf&#n^755r~LBE+WGgzu%MYc}9_# z@T)OvAacw9eOaz9{CBX@FZF*6Oed1AQ>QoVU1j{(RSfksX|ae`f52@fb)^4c{+v*V z@TgScW=4OJgj@t_O0*hqwEg;0ECojq`E;b*_Q%T!J=^wD4-xVJgQ?OeJ4U~*so7E# zVW$icjv!V84n_Kd2E9@<1TIChvaoN7KF3Gu#T}DuqeMpO3Qm&N5zXZV8P7={I>d2JZA)M=kxDcq*Ds-e^MGmqL|{}A ziOK`c$XzS;*|=mI5<7Ygw_j8Iz&Aj0WN!!OgJ%+bn>4;o)cA=qbaXQRRpYI^6 zFG}BgZpsTx@h;tXF!A^B#`o#A2O%ZnKKt{{nhPeQM@HXlYqlqRSsY0?4X(Bt`B%-% z5fl;5;t*6lNC6x$6lj`iJ02f1!Ac_-Pr#a4t}E|oLH@smMjF^<@tSQSYeG#;x-)}V z*|+0-HRsaXtJ`vDPqFBkGiL~k9C`tQekHLi+UoA)&4R~M3K=7GA`yX_#&{-yj~cy3 zry13FWUC~j3|HH$ODW5XBLP-;)i)u6wcc^8vQkov=!7(YP~3uE3QEGS$#FB+ObgQoG;p4R7e`?cRz4mnyG{d_bePp&o2EnY^rmWamj9$@-KPK}p-bEGeLZ8p@QD^7sBU?$pM1WFb zD9^)d*GCxG2ktQuoHDM#9h3qoL@~Su1wGP#lNj5x2F7C5*>eUxJ@JGtUmvmQC798` zs?*~?H+u?f@uBVLNRhpMbMNIq7d%cpNZd#ui9M)~(1gw#_&{Du`tj}%^Jc&pOw`j5 z-C*iW5VTi!5LGu0}}F=bkkmuCv@|r8Lu~6w8@AWf{ixsFo@jPY zl~K@bn*N;b^Y%+YioL4DZJ=&GncLIPt3E8 zt7Tt!gAfe$S)ix^bl+9SlbsAWbm6Y5+={s==Tr{-n=N9l(ZRXZ9M-c`pRBR2a z&}m(_X@BqSfGz1LB4-;{?8+-2Hzaa0iZSeNWS2Z}f1mkZBmG^cI-hf*cUg;}FTwaY z6K!JROqc;ecsSgDe-I<6N9%+4knr2`6A$lNCr0Ml;fdF%`ZLxbU^#t-fJR7B3k7WK^qu?6 zkT2E+q73sk3aee~^%wP9JDy?(O(9owRE|r3w($Lce)0GJy)&CR8e@ zr`w|lPwKCbtp=h;w(<T?`wKfuvh&Hd_T^1JjNF*#GbgWwt~j zL?uJ6e9c=oDxV3J9bYH0_ZcPxe*J?9fdF_#{^Xn?DozCYkkUK8i4b@|`@Hwv4H&SX zV_z|?-!<~9;DwilLpK+~1k!Ir3XiEEj$d|uKzz09->*!6&2t7EGv)Xl5yQ;a>@kD&~G}Kv7fPCEiLp^o;L46-6_v#Wb^$)3$`XxMtb!#ybFqm&vn0 z8lrZboUJ0~h|$6fE(|cId|BTf-A1Dsd3;s&-B!EaAnJ^aiXTQ7e;GYlB@|K}NW4VV z!L;}|`r#TbC=>a5=#YvoPV;|n<0B&wkZ^Tevw=MhZeG~%VI=#XQ9Y~marz~4up{}Q zqVzYVG(hBvQ6AC0!}^5|#yFK?=I67n?rRu8ibD=JY_bXY(GPm#n#~ATi!R>=LLt)F z6uanXLatb*O7g5DJ0W>n|G_ewQ>VZ*fofz*z#6^WfZd_=p(YS)nuI~jiT z3Arc4xPfwud#C4TOe~s5&^kt#oelfddum@Fc=hnW2J+-0UNFHG_n#5=Oe>;lYFJnu zJ81UurujUD()Vv^r*HGh9HpGUWBJ3w!wI~u6w>YQ!}-)3vx*?`z)gulgwKG`u0N+w z+%|5%;%5htC`sI8x7U}OmR2*;C*I!Vwk!q0XevELN)=JH-Wt)`jf=9NRwBBUTwrWp zHQcMspy+c-RPxGxAUJ9v%kZ=2j?@9h<(Ep19C-D!hhPH?+n8tdQsc32uB+kmf<_wqW7!A zjaA|8*E?sDU6nkB>HWXXsaoV&MH8Lewxe&lM|e2|==%kfBRWKuY#TB`eAJBB3zhqE zsRLOL;?xH2)tA=sEX0{#Yj^Af|*W$oMdN1?zeAL^l)}bjd#}Y;B7GF4Ia&(tv za}DY$TVK^)>j+^GWJBnzLvRKY1#f0`aH4)E?(r;oEOU;!CaDT#`#J(+50pt(K5mwe z`i6}`TkU?bteu*8M=y_8wNm@Cx-Irt!(>?eftw89g1LUcwr0rQwwk7+qHkY;GsqTc zOkDy-w=GAuj#le|_(tCAaBt}0$c1hK)e~qsUl(0ExcXc8r_Y`Zs!rmWhvpa-w#-h9 zj(DGaN}KZTcJ;ss3vOkmq~9vetB%IIM8QH!kQjx>o8>DtvJud>mtT0?<<}bRnk`m# zPDu^D{uo3#4wxxXI`MAjoU1xu10Xmn>wWn%UCMHz>D4P^-iO0;Z9ztI2uyZo%Ky*^ z_d6dp-Yy2q2_Db~OgsB&=Y7d_C%)2tClzY6KR9+<_vt_Q0N7Xt4H*I!R-xY-#5|v7 zD{~INp%zdq7qqxRM$4eO|7jT-dH^_3#+dvXzxQx_TgpZQJMhy`><|Mo(*Boijk{ld zQ0em61F8@~N@>Pvzdb&GJlB?W%&hSTZ43_B9s@+4;GFIQoF)0Epz8T?Wl8^)sSr!1*yD6(ReKK7R4VV2tqK zSJA@I8>3nnL)V~C_NtcsW5Hz_NYbHT>*5?u3&5y4U}0PounXWX?=hidaq~7i}pO|Dc9>^^=(=`8WJ%1 zulKWAllr(@cji(WM=b68sA<3b)MVpO`Lta&I>RnL3gZC5RXF_^fWZ5QUOOJ#Mavk- zQ=Kw(8eJ!}%J79_FU;q;;SJzzd{wQ0a=2$HW^j@O9t1z19&S;}RWZ$bT0lNpcFTv2 zGrOX()6lN?0~Uxd-}ZXj%lR8DbbK}<^weO!WA;{F5y2Q|X`V@Rdrj-Lb{&D88%a88(w#_%&dA<%_|9T0D(>v`G}mQc-;!-# z>?w0N9Zcny3JZzc$^Sc+cNFd$bx_&68CBl@g`_I(p#?KPShsD*q{ zl>@`!U>3=2k+A&R&dD|9uJnw9d1GO=R8fnEo_sF7&nd z1b1VSO#D%qujA%qRLfD-?Yzhd9juz0&eC2Ra(hlY&%#%QPH{YT{<*d>0BYhPwdY_6 z=J_-W*&(hd8a<3WG@5O!&X(hIeb+=`n~b}8^QQ2?Ub+KjEO3+>#8e;R|rPbWKX#SnP?5CQ5Riu&@mrrNb#Wt)Q{qJCJ*o%AEqlLS> zIA6h>MJ(VJcWdVu*=eE}sJ$zB4SsP_wl0+6(gYMAW!_X%V zHFUl(rnaLaRQQ$jpg}4ns(6`$yT9ubJ?rS=(~YHbpQ?s_PgX2$!SEO64WIJ0%7&h& zyG}J(e`-w058VIVZ=I*>=iM0EMr!cm=b2-qn@A@Btq}8*BM%z-g7AelO{qs8_1bhJ zy4R3@+7vTDVJnZI>4swouDby0_4TOd#zPW+Js_K8`4+$UfNMbNp&H5rFG8maw*vV&a1_}yia zJN^&sOr|g34>b@?jS=Mo+@uY-xbkk_cTqHO>XSd*+$_{o)OBJRAuyqSBho#5;VlYs zvI->%Ur-E#%Pnc-b1%&6j4!@A_rz_PTKxI%a&J^cSX+&A$ya>OSm~4F?E|5H{vSA< zV?{(`vme)kI#DT=k(;o$3?Vb|Hc zA+I`i$2mO&zDDKQdCz-TrLS<-l$YyfMtd&70nIdRE}_;pf1zu^8YFo1B&F>2|S5^9n)U8$npX9g>O--Ecr^r^^8#_k$h2k5&S-^BVP~q z0=?Kor*t?q3@r(=(gwGmq?rv@mAUA|PY4O9t{@e`!=L&FLm7SDYq#;n`rGyT*gbT_ zpe5hxiu}x)0WR5JL|JF07*5n(xxiKJ?{o09Rp?U$)38pVU#TXLNbiWNV-*xc;hyHx$^Ynjm+e^QBucduonP^*sj zc&^a9IOVD2R(JghU4^;&jTaAndLQrF7HpYh*7zL9N|-I>!X?$QBeWT+w`t}-(G5dz0yHka;?uhN(*1SCg_B-os=dQk z)9r(fNZ$?-Wl=mL%YQpvXE&yJ>lk>QTCNj5-Wk?zGtpDlxR87V`N86+Ee3Z5=_weW zL0qHIg!C*ta7s|KjW##xs@icvU8i1sz+>pqVlw_GjNIcb7K{|#^eHl)Mt$ae-`MQi zdB^H~g3_ssC*BCXl9MZQ3p)L@aIv*y)tV>qU<(d8O$FWBc8H8w8w;(Csp+|++bW5g zun)nU?Y$UBMN#~t9G>^qu}tZ(%GvU*nEI(0KGqXu^`!d(o-E>P#16c5VN;TW%@{u& zjh);iYC@!4q=|09`@Wwy-e(zcE1demfgT$z1)ZkVkYwH~9Fn}WeZ1pd#oNJDjpH-g z)ht0t=0nt&7vy}H7IP>X4#&Wb%ER;@Xl~6LMCQCcN9`D#QjN^{hKq~l5`En>mi82V z{x#-JU9nqnB(ZCDU{6v{Vd!1VSHv@jNN8cF=;>w;Sjs}X!gQ)JopyF1Gs~sx_jp0g zdrqVuze{pX^m!?1U>NRM3JlG^r_v+z~P145p$0s}Fr3~0HbhK}6gtqcN+b!H+ZcnouW3Pajhnn5f{?yDP zwz}0uOj5@*LqGlwg^H%dWpZA+Fg<%%ZFt0i+LoUbQdZhm2gHLe!ymbuOBVi*r&94tu5h6x_K+D!%RAFWsHX}fHca@k)0 z)X&~R~zG<;tCNsmorv7QYxO#F66(eh_ zs`&8z*C~`x7KUJlpcs{~%6}d`7l(WKB*UmPCzMLn*PlKK)`6wNe%;1b=oh{Nd{ghwFp40qRzc zYa2={c=Z>~RUkC~lYPxYp^G3wrB-IHJ%UMI*Z>)`2~GT=vi`2> z9Y5aYUf(1-$G&M)GZEvX0b%>ObIF7H)Go|9n11V!!Jj!<<0IcOwn27SSo3-2WJl=T z%Vq|*RdbK&Vb2c5o-H*If=QP&zEK>5iY}puZd&H2(^m#AN|YQ!Aphd@hUouIcAHmK zll}7h@hX!Yd9AV1o+0sWfmwr`{q7%+VcMtIN1hY8&kyA#jM!bb>YE`m!?60+6=@1b zXB2%m-)egQ#g@T#^n0VbPPLS@dziE@%NkRXvL*bydM`V7$hKg-%32=jEA50CuKyub z$w`gm5q*10s2&PN(yC%cVoguT93_yv8NBbliTJDke6EapG*y*Rx)Pp$sLwl@6Ywx# zV8)%Eno66S z6hslO`&LhVPe?RwKCUe}$-suEQ^+9O@BV4{^?Cl}4TU$W+eG7hE+B0$G` zQYfSaM&4cGHvx)9M}(-lVyk|qgn7R96+qU5$%y)7emZtl*wRx4IQ|5_JF>O<_qh*Y zYuMpMXG~+ftK>`bkyNRqsydrfvI31?s$(>F$xQ5!7IZK3x>4 zxPDhqQC9zj#ih>E1Ls%yTH5awzrs$VS$yYc!7W})dd&a4q4-g0pwOpS0WEz;B>I{= z*wqPk-IREkxEVu?j*&R{l|WJs=kLh9U=k(De|9~+f%C}<dF&^SdXXxl5uw4(}F;(N2!rT9(wmCF^B>=LPp1sun4JFCiD!m zsNWL38_`U?YAP!)IPHO4Vx@){f>DU|pDym==jGZ>Z(xYQ+qQ~=Y1scUjUbE2k0JW6 zhYXx_t-pDGwj}}%qTnOy7Ne`unwWt{s?9_Nw9v_TS#%8&ca-6<(;2+03Qj+_yXFGD zGP5{6Ov{!9I4`Y0%@1R9P|={gH4Q$t1Rao zPAr2e34$f+OOOV*@4a(Qt6ETKWz-%%qZq+c5H%cNUvm%d{4Jl@M>>K;rqf_K$3_Fj zCoN>m_)@&_=sk5ddcPz^UueYT7TxmSRS$&(vR8e1j>Rq)bu9=QV(qMPUEUX;uq4t! z)Ql!`oGSheO2PI%UKh%0>R2S5(&%$Ydd_I&WYN+$M9E2j%!*(FX4W2W1@AqOu-wck zq^AM%#cf&fin#+%Bt(S7){PKM!j`vt!Dhw9jD~uTy)kLcLNUsRq^@7VLor4h$os`% z^uqr0kHRr;*$pWm4wgn6t)q{hneoa-iQB|$s216lx1~{$u24m1q2rUsW15=I`jtnPD%!z)BE4ID)5d$o7+RCtZ-7gM|)~~ z*@)l;x-el>VM${Zt82~{2=V0a0;nX=(HykICl?l@r-rZf+rXDXXUYc*_5qsOREtz3 zh{g({j-H%?fMx$f7ylHypc@9R2nz*;>urP^ht4VkWfb}#At|YtvH?h`m{Idtfoe&; z=7vk4@`MdGnZL|RTBfmej_J)NlK<^4jYKAEe1G4@J*ivk@>fgj+_~YAO|cdxl%ao! zJX3e<+gpvul03{Futc4M_tf6wT`qqEJdQy#5i;OBJ7vUXF z>@i?Jg<#qt?N%70i7ikBT{W!BU^K-s5ew-)x*N;{UJyX9TM%A>scJleSzE~roQ$HC zc3aQ+8?5Xq)2Y!j&qAH6_nco3v}Y4HyO!xk@)M!`mN+X)Zy<{4E#g-sA8Nk?ASMVw zHx-5KJ!lgMiHNMi@9Uu8xqRiyfg}4_eS;FJ?IoVf-jQxXEg1AW3I<=(*+ZH}Ow~j~ zf~9u{5-VtY%DWbMFcSp1j#nAVa295bCwVuglw#wAKQX%_#R_{E9C9CjOP`AqD%3+h zgO6^3cE65EEi|1`8%~!}0m)H+M>FL8i;L%_t&8rqgVxh=GEC9`)_d$}B;ZCR(*d%( zU%=a7LSZ0Yy_<~EgZC+vw74+eV*29nk#*9>GO=$g@#wN+%8EWOVw=QQjx8b_%tg|A ztyRltLzC5i60LDkH$54=?3>^&+ML{eU4|YT*{19&1`u6wA zc1j;50Lu>4Ck{l|i(0>3n!+COD|X#Z$mIst2&GHKPL(2^`J`=h?7x919*!kn7q{k) zS+^&y8RE?NE~q7pvA)bwGP?i!d$5OKy8>}`t8`dpnzU8EanTeVDVM@it2?m8_72;V z(&uyAOZ~bv$|1Q~$O}#2Sb3UR**9S@=dWbDtKnr&2tf>wiNh%vsrQ1lV|Li)Af~I% zsvC`g-qKC<_#4l!t(z-qbD#Mo>HYUD%h->T8@_yh&m!RePgx;25O%tRN+}^xU80r| z3avG{6L>$HB6*95@=yvD@1F*0X(Ax57xMS_8&221u57YzxjrKQqoUBIvv3;%VVb>H zavXkbW4&^8QDC>~9N*vhkwv1M1A9^QDD)mlVV7%HuRJvRe6$Qc)U)XV5l7WjHXm`P}fgjX>@Y%hl{= zoFI3|K|&uDf(kHMe}wagYSks6X8H&h02|zjZXAi}EX*UbWtZ7E&t_6ATOT!Zuc9Cs zyhO%D_g3VR(&h}BR&;N!ak67;0->hF>+~gSa-{lEk2|IQ}eh+2%Xw@o~Zd9^g zLLcEV)~ky#<}C0KR}uaeFdrG%ImnJT8)2YwvatQbVFAgJXu<9{)IucJD(IOH9?OXLDrYmElGNJUQxO6_;4yVZmJ0{9CBsbuia8s(r z8l@|rtL|L!krNAll5diq7<-n$^|#@9*%TCxOjYUM?n|J~gpGNjI}e2fmLO2%_TKhKPd#c}Drlw;V<;`!7y2 z4o#nj76u8P?$gG+5R`Wei$Vu(`_gQR1jW}r_M|d2IpIoiX}o+fCPw$@o)JIN>->i6y};mjR$O$ccf7~y+Ega(kkNO=@;qr^fXqG| zq$F}6n!2?7c@|fT2cVm~W*m#>|MlsQ4yB`>}*{7CYkFtaSC8NEah zCx7G>>mgtQ@x9qFCs6J{8bHzpG^tjn?Q7^tPDrX&yYV}2pfBYHVhmxvmbQMNn~{Rzj{iDJ-L zQbFrSTufbI3D=+qz-%ac#S$x2TaI^aprQQ#0SqY{cAjU<&)BE05H_Awfnf0))EK^n z#1pgTeD3wPKecMgtwtKiTuP$7B|422pScW4-S`qALFKt3zL&_LJRm(Ibl-(7b-fq6 z?)WSswp1k+D9Ca?<*QlaEOZ2F-+%oy5!Ml|a@}@M!!mxM;R7H35l9{W@hzOc+QV5Wg7c%J@h@s1 zR3r$y2s`;$o7g9acF7GVvUI>08!c0wW4oEq^*?j=th{@-28Psz=L@SLL>0U4v|D%w z*v^2;skwHD_WeQpLZS{-Wiw)=>KsH!|sh}A?p#)B` z3yj5CJc0{@x+c=+1Sp8^9GLf70}VhpA6NnDT44-C!=rNk@0$^eyY~5jcl1C4Y8d|> zxGl8E;TtwQ07>#F;nmlXw>Z}B6rT=VO`z{oVKQ`dIKVBf&J&qF;Z)HSjP`SfUWQ3~yk!WT~V zSB6Cy*W4{h=^^E?{H`aag)hZr<@Ik7G0J(5+WVIXkCC@28O4NmbH(4~8-hl{?3Hvz zjxGPXS5Ln*Q8ht-`m+ENWTl7>UBS)6bm0}y5FxRA{JrOuV5Gy98=&^v#qYyfamVQ{ ze9@GNPFn{hD>4K(+cYcs%sj{g{!{7fE5{coBk4uY*tFQcLZey$o2&Cd=iMhGsBU3u?)NExe4=Z~{@25@@I}&?9MQJTDuwbc3#rRp|9a>Kyla&Vx zU;d7D98!j(wBlicdkBVJLY`}{zzZ`f34)#(g(lQqxqR8g<*)VD{{^OrK(28@oEiA` zXkU)u5ojKP2f&D9_^4evin6-YTDacrdej@tgEA+yNhM$Lvv3}vzTdT&v3QFUqb#5e z%i6TL-N&NvC>dWV>w4PyZu|VZf&s=I?qVx*2eEJJv>=Tk3a^#gLG&v^XMq!;73|Q_ zN8D$T{ff3_+(ttNrHNoc07CehXz#@ zh9i_Jr6iwqj0Ew5hRdHJv3ZR=t{lT#%V^zZeG?pR7$f;prtr3OyYbA?2P2tVu80~u zi_bPpc?T6yQ_*3t*zmm-3lNwN&hZr<_G9yVG1w2%#MN91~R}Y@D zcV1vsu)ezJ6CAz4k^9u3pf+qV#vsQSju5>WmNQ=Xuz=Kr_bx;G{|2t27LZ(N@IpX< zspF(#WguO}5x$dtuu`LJgKQrpu<`E?W`ahPtz-cm0VgZ*G97jX zT8L8-mDxk0>+5NZ`L!>}eoIgqh*)IcXvu9g&$W{8;YyW(q!gPvhGygg3BhfsRz2-p z`S|-~s|3HZ$+JQiFp^m_6dDtE47N$?vbL|`*MlC$UQR-jsq%JH|M*4z>XNv16qJ!G zRl=vt##B{Zb3aRb*Ax|c5Qz7)RQcj!{Mn51nbk|y&m@GOySyljNtHN!!lRfobaU#f z=%=o>QrwS{4b~PHXWj2V^W~Oy%xg-_jx@o>*zc(;a7Qh;>+;1l4If>!xL!(^*m0tKN{4R zZm$HL1c>(W_7cbMbG)bv9De6T;x_-8vaX&)qb=3WDYoH+(BkRose~B5P7$WeDzTJ9 zANoT?SsyZ|A3gPmv%bw1BrnRlj z0XkUa-?j#){$7=}(%2wUVei~o7p>9gS)&!um5P#3c7Mrylb_V=DJQ1*R7&r0n5+O4&I%>we3Yk5W@272ci@E`g57<$Cz@=PG8FmYM14R2vjeYN({u zLr~nUwD#e{hvbKow4loRSQ@SE=YuEjpW$1LC4d5A=)L23E4w{Cw}1Gpz{|^9d!k0J z0HpM;Q-g!(+qW*MlOJYoo=jix%d?_+{g4U@)%8H#MI$^iwl(Jhi7ggZwSI*@CK_rb z$cE28mD*vRH}u?9*y)mdeY6{LGOJ=&PSXb*hd$R~D$cK+$~{66e8F(r*yApiPJS>d zbIu-y$E1g_Pf#6fGglSoTU+ww++K{1-MsK)lR9yL34L%qaRS(^q?>@zBCpesUzIUPLuplD?||APQNx3tF-X6 z7X)7-91J}X>I*`I@w6~qQd#E~Fn^HLze)B^z?q&Y(NM%V9%h$anQZjW;wN; zClMupE06D)O3HtGa{PMldG%_0*3><8zv}M@@i0a2dYvEh`su0A8i%{456UOj754UN zWGi@lR3h&%xg5S$N4RfjS8Hx9KC1WGQ%Cr6a4=JYhfYtv)PrgZROKWB{%Br#di5YB zX^%LryY7+^XWQ1I@#-K}a@)tx*p^6hC+7`i2%fU$VydB(zdRwwqq8+tM}gOT`d_ZQ z)c#!bH65jU{M4Dx8hO zX=wGMR zRk7<$=H}*6ym#?|PAbM=|M{q=|3~)7SH+?Z>BQHfTT_`5U*Dzo$!L+l<8bn8`0?u* zKMg)U;es5pFY-_N9>JVx?u-e&Op{P8Ast~hkDKV|A8uNBw`OD)M|<+BhaQ<1jUToi zy>b1X5c~h5>pj4+Zo~KSCnK_woh?#UQe@Apl)Xhl%FK@JojnrOWA9C7Wvh^6Wjsap zCYx;j*T?(5-{0?h{Qv*s@H)KV@jUl)-`9Oz=XqY|720&225_F8fF8@ad%2^moTg(3;Mvo+Npn*LnCtse zfB9J*DSL?}N;*&H1zN!p9jtbBnvSS__-8fo=j&-pHF*~@Fk0BtLMo%|jzQOz$=BTs z`6uG66O^>?3&>d6=N&)D^AxZS+E_A^uabd|{t4Q2vM>MaA-yZInd*Ee(X87C>PHKz zKHv6`(U!QBKWB0r;uT2(=}U#vO=uKz`!Y@_ao?@&+gv6)$A`k*e!W2`8-AyRG6NOj z|CCXp_Pw77& z^VxG!S$^f)VvGvxah7~GqFmzt?6L?Ktbt>pQPS=`bYjrwiKOh9A?7v*PwUs5=Wr!# zM|k`A<5>!E-^enQ4j<<28*mJ(WBy-R#=gH5zaXEq=Za%s5KIQPWRah?9&(m-)WhUI(>s8 z_;eip2?%gI;RreagoNE~fx8iHZlBIg7sG8K%BjJl>vAtQiG=k|@aU%-q}O;`5{ z7ImA$Lo|CqVWEjW65RE zk_Yv)kd4W}>9+KfLV^0$wCmYS@P99I$)s4ouaBSSNfaBp?)6_jBJ6ZEeAv*zCQ=xe zH#z5wFc?xQgikMO6W<%p7x z%%o+aRqt}3@9pn&puZrbzJ^^o3o#KVK>AJ2+lu6x)ErCGL?z9cnv%lP9067{%nW=7 zf7ZJ4&UPgimRR}uQ&hc8v>g^!9-nTKJ02;uq_6LKUtE`*3XoKKS=ogNdm7Yu&9!SE z#<4(MuQO6vWdF~RrZKzhZ{^hA6!>HM;n!8yw}&3q!vzQ01pmEGEqC zZA3hca-RXQ=?1Y|T8`_jn8nWjIcyPZad(8Gu;Ik_3DfF#m2)HUYRC|_^?|ykDD76- zkN-J)Zxz|&AX4?XFa(=0(?(*UxI6qymO@%gP4oUJvR(QcE}z-Mul%R&e;(WB>n2Wb zk_^SJUIOlLwHG$S^J8XVse(uUt)|XdSSO;-UoYOdMiwg7$B%<`20kp%vvl)nX*i{! z^i%&vhq9Mu!?*aA$-R8?EE97q;ZKIhXF?!`J8jJ>EYC-cpO#WUO3V5Y%~^isj)j|A z3k8N4mjPKM4}d7+@ExqfJ$6=FfS4Q5mSe?}^KNuS98X}Wub$Qa|6b^_8&8I-^N7oB zYb~VJ>{|cZkmlFol~`gGnPU~#eoNy{{=dUlgf{(u`_qi_|0B|ssW{RrD4Sl6`u7wM zPSLkrhzkUldUz*X;;1ajdrLiMxbg}LG<01HVI%bMBX{(_YgWdXYe`Gr=4V0!i|~${ zu-djTaRyc*RuKpAHX(jxO$yCxl4g*?{3$xEaJ!6f6Z# z*N6_WmI71|PF=!yPm&zY;>lbN#%38TRAc_4I;@}XK3?pCtf3wwJX%c&D=9t%V)${0 z4$dMxSZoAt8$gE)wwdx-pxLpmE_%oc$YsSy+$j!~&$N1(N4eOSdme)y&rA}iS1RMz zQ?U6s{doT7z(B-yp6iZJ2y5zx6yy!mmZ~`fQuv%B9L*&q?3-5 z<4^`bLXX8<;DUyu56?0L@9x9XyX$hZLlRX0Q4cdgwJj75k${-PTTmbJf&Cq2|jU5DuFXtc$MOijqDs`@+4I^SCl~KYey>c+3$zyMPiUf)bpfMxJvl!3j?E7WDY)=7`%KX65PZf5V$^)4N*Asy zo1xh=dh>p!`d3=3Reqy#}<9@zw<(|iB) z=U9F!=T>(UB*I>IqtCItM2a6GBWnAE=IoD!O8l7V<=g3MT_PgUIhOcAGRw5R(D|K% ziP?ZhAcZR|QeLh}daj?J_?Hb;4h8DRF+dk1Rpz}G6p)5M{@_dU;2}8)ym4oVK1;!B zU#0j2Jn-1$VZU?itoxcZ!B-?G4(cOsQTm`rUcZbFGwIhoqdxxmnkYT)U@gvg^1Lsm ze2}{H3Ezj*5zbVTN23P($hK zgg51gpONQtZrg*N#_-La`oTI{fxXc}WeS!N(QTpA6Vh_J1<)9ZdpdubLI~DFWiLHh zAjrY!=3h~lWK{0kn3Y4yybExLy(#`a-+AV=qNpLnlPSu8fBrJ%p#Ctof$K&={QSK! z3L}iuep?8LA<5D5o6pKzW0eJ8i)vsde&fwj%tGuyTq^evq9i&zv`0;$qX^r+S zTiA7;zF#-ES~a&%xDlev&+6oxzF6r_>T#mqYI7`kHl#cI4>wZou{&X18%ih%-zVK0 zxYgT`U>+)_wkA_0@k-Lo?|hwvdix_X2w-q#APv3rvmGqm=U%TSo9j0y&!1j)B8A^v zjC7SB6&V_BefQ*Fal2U9l^C?2+5KJrJ#>F3jp^?oyaNh>El|Fo0GZ55eEP}u{YZ|@ zf*r6HRvP;5t|u>qHbR&9K8#8Ishhs$A^_pwfz|(=3~~|zc5&OxEYm`6>VIc^%wOE& z2eHxt37I+Us_BmQ-eM(UlAs56f&q03--z%04gla=d!!#DM{ShF9)B6RE@#FO>&da> zzk&!GEi@f-+kP=Wck=s9r<2_BzY|DV|IZ_{3_8K}XRZ&|0u+#>_*#$R1d#{GK6Uz2 z790u0WAY@OOL8KCw7v{PB9Crd1lkjYkYUjAUAodr(B!9P2^C31*|pr~n*_ZN&&AW= zz{S9E03=p4pkJ9aHMG%U8~CT_2m;&yyS83Cc0_`(R)D))GXx)k#gd44o{|fu^#mCF zGL2RNcsPwyLm&_soU+rc7oX%6Iw79Zz(;mnWDBXXMOF~nQDP6>qiGL9sdL9SMG|!zD-B=5^ z{cv?`=i>>eYR^TA23X-)6rgCm$F37W7&yHYBlWo1gFA992@cscs%`0D*Ny1QI@CR% zfwCWRu|C}0nbL}TeGp*Cb_%Z?9nML%-Ru?%U^*PfTP&!1@O*|6q6i#0PjEQZ)vdq+ zD=yk<+qt8ND{*)D#N~9~+0qRR3_wMusY{vl9!KUM6k^FVR#EnQ@4yJ;&v3^ReI^ZZ zAFdv@8GrOm@fXW^_m73M&_*Wo%k}~`e!byv zhPSZ&6V!R{Kn|8v+3j<@DFlPMOoneECK5?5MPdjNS#68RL>3FeKAom_gUV!~vR#2K zIhJ9*f$Xda0q6Vp;XaaUe)|BP^9B37g&QJNkw34hzPU$O{5tv$x#<0b%TafRcav{- z-8{NxryzaSsL*!!I??T-+@&qTE}LPh=Xc$kFU}ZA;2o{F8WrPeY^+Zo$O%|ff??Q9 znfNrSjG_*FD1>HSY4$%tHIb`|`2qez2Gpgd?!+rim1dxY9YE$$P^jmJcmM&UwSNsx zzR-MHfI$|unLgCINR0|YtqI@-Y|0v9vW*+P7iV5weoEWsP6?lsE(%5#ifQ?h1j#_Y zn6E{E0=y*mb*g+Dj;|z%H;410Gv3` zL-+)o4W_pw`QO6D*W!}1bha6;ma7nWeqt{a0q6-9jYw);sV51NjWsL{*yf$1t)qRZ zcD4V~$v=rC%JD(sgNhq_4=&f|xQ5DRv#hWs7FZR`V8M;q`W4Lf=y zXU`eE4jb*(zz>1B950YX!DE;{u{i{U8ig`+Y`-kL;QHGE4@!29+$8^#z0{Lt=rv8L zNU3@XY^aq+)G9CnLJ$XP9(~y0D*1qASHDe#8iyXtwv)r|6XXkBL9VYK=TZvU<0(&6 z#gQB(Ox|`e+Y*g`?`n^IU4s}J^CH;l^0dKu&aav*3{dWn4fi$xTh8u(rV)Jst`GEb z{sA`WQQ}qJ1!seiEah}H#r~YHJD;kTlh~bY8SWFypqPPSFw?n}|F&xD54%ZCJG&c~jeHFyS1<{JJqoVEJM%;7enWS+6(L|3ot&CR({%t zC=P&OBun+)&#wuigi`verGx?ZMgZl@Quzup`+vw1mPV{#V_Hi2yfhi*3HtnK0^A!l zf@FX$oR&3r`HkbEc`#V*eku@Ms@VD;DaY-v^{hG@?s--^Qw&^##RteHt)E0#!Yg{- zJaNiSZp2>5o<97JvYS6Ja#bhl}K#K@|KM8^kIrm`xjPf#Mo_xAKKW5t8@G*nb(;1bc}Vq}BRAZ5u`V2$_nafRLo`}7(mipBg? z(wV}N5)yI+Q|2;sl38+sk6MdY zeqkYER3+Pau|8*Mke?XH%_i`$idD%)I=)$0xr2_W*)i%8N2 zkW3_2S6e$BUJp9w6id5(8yiuwf&t9_=#k`|%2d-iDKjFIQ z61(}{ZW)@**Fh9#3j7J8r?dZY=ZXmgw~)3eq0R3%9e$d;n1DHU&4NH{BpvBlmUG|r z%o#9;wct4};W~3Zp~;<+k=rl+27D?jVL5MuD$P$dw*8LiJbJIhZ*Fz$#CHa=p#u1;c5DDEQkxpoYWD+z%-c zVNW6#WP}Sb;5YgS^R-hL)q>2VTkJkOF8%!)E18#1RFLU%cHWzb28{DHC+`J;k7`HqBO8Qo>jFN51! z;FOT*N*cm=AR!=1WbtYeI6n#fdH4CtU&(`C-S&`2C@MT`AtQN_?!hYZ40mx;XTgsO zmgY2#znfpB7d17_W60#S)r&}gzzNb1Z0o6B8$G+ou96h8zwe&p_4gO@fXEhvb#Nlp zn&s}W;v@3B{4adL@Y0TgQ%5B4rdTAwWBr6xDXM9h%h=j@a?l73gLbXZSK)wICBOSj z^7(=Sy!FEa_cPP6FX+S$XuzubcvJF(96G>xeqoWBZE1Z+kjCCR z)eknq@}ERDJg2YQoA?%DB1aEU9pLf~XQ)jHqq1 zXV`HZHMtv6e(1QFyPB8|+;BIsYt`!10!Fh!Dk>iKl3RhQ_1@b<%vBe1o%Sfb=SBLh zgPgy6?k$G^6|f>m`E}Z<&9uzlYMsW6Ax&(D8V;Vq#3g*hf)wE!P|zKMw+2Dx zRHT3uJ^<(y8Ga;-Xg_Sw;?LH!xDS-wRHx}R!Hx`Go~n=fuzPHa1DYR(yYvTg(UV|5 zXku}_H&FL%B*EMS>Z^==MoNkK-&Q7H%8UF{j@zXpfi(wAY-(Wo+g`SZ>~l@IOm&Ck z4e_+49J3!sYMXa=SjHw5yvE>;fXNYe;sSJ`Ozq6G>_=vWLa!nnsJ0Qu{L5B-mf}@m zF@1WN3ftPdb9)Y+j6FJR@5!s5jx2g_7aF^q`4Fr4<>+ASiG?!qqdE~ zHR%8(lEr&fO6DU13V{3UnGTJaAvgUV$2K2E=w_3?j~XgiXc8 z#V*~#bM%4~1Ro*6DLLpcKK{WV=Jw753?gfQnOV7WC3%{i>>S3$Xa{mU`xnG?vx9?! zH8mx*pci=74fpC!{my8d<9GMk{VKtaZo9Ou1H-)eg#h##p`LQaf4c+l;L`?4&<2bs zvQ5H636h0=06@%T!~5;M*3L&aOquo#*~vnaye8K_6k4Kp8GL^xY*0qG!H-|6?&asX z(hUlUxi|1DP92wGL+GO01|I45@e49Whu7t|z1%^5#K&=TX6)wMh@V~BQRLBWPEqHY z1REil9gvF$WR9K-@~q#Y7|{~+pysocgQk|xm=6I=CkZuSe68+ZkeiiChqnKutz~xS zq@2v9Kz<%cNTYoY%@Cyo1ep)e>2TfW*b`~I0+|$(_jM??d zu#r0rndAggG>2~-xevjX#0_k_k*Qorg^KcJI+|B-gSqR-C~~9|L9f!{0+jdxFs=GR z<6fuT8q!#VJOU``>>(SNI@%d~VXm*3t$|2X_`rBjprgn?zY34dxaQubNJx%GKGnjO zX7L2d*I6VvZUtkbvJ$Qp0x=rRrDhv7_K~4YnweBB?o2+_)mX)8!#E&|o^OekFadNK zDm(7hxvu31qkj%H4z@kRl}cJJIfI?jDxX==kw#&P)NG{g#LCgTHmX74ql2ApeVRhH zWP-qP%0Ex*w{^hAo&gc7Kzz}JR#wVHLAzGIHu`|$ql3(X0s=7MW{VP@rJ;BF&bJGp z3yh6ln<0Z*!FtsGJ6{{(P}~N|hy#565HtZOU+4&#cw-jT4Ey!maP>MAh~#aI9^Bj= zpFq?f@tH)F2nhnntb@vQp@lJ7B;%#2<()d;*ECC3beVvIet& zF%h&buoQ+uD_lRfevGXxbmcJW>tWKplOwhtCzAxC#>+ddeFpdQwc)X2D#T4YGg=}H zY9fu1dq{8G@s7_4H9Q)mg1hh&XHn0s?e}4NIMVNkdIz0^NCbz?C3L(y#Rt+&da%3M zojC6k0xy*w#uT?8!3siO?`jC2DUEzMZ8*f-V2bT2U6caXLOIlM8Lb&rY2jZFEF`?% zK#px=N3-VH>^s|MoI~25+G?qAU_u+iOcgGT*-kow>+=g;cJxr9aWIBX#7M&@=4J?BJsmyAz16F}KR`y|!t zR>N8l??VuX?+%oD2l4vsb+|#r;lac25vuW@ z)R&_?jV+jQgvVZ~78qstu|fzD(!G@^C&ozabR+1q$(WNX&J1SjY_VG6wNJH=Pf4Qq znupgz3JN&2w6$3lIyn>@=jixDTwV@Y=7}ScNowC(5V*!xgOfG2uLoGw0(Dpe3lK@z zgk&=#pqcpqYHSzq@qDy18P5C62;(AXH~Qy+=X%RE)DE~CBLA`qt3~a50I&PIixI}w z%E}Qe-c@!$KkUL=`{xESu?hJxKth&||2&a&f`RCZ51cv$R=^{0rslzwMLx6E53uW= zA>WvZ(DtuBdlzWyoR%hBqjo!c->#QLy=r$u%!`+4x zCx*-CH?-PCv{_}5>&OIDb-SZG4KMyPfPjuv)VQ_r0WBi0gIg9cbKCjtbY4YlT~-XF z9~KrCQeIx-aH)@X-UYuU2_!BE!Gr*0aB&38NkMFmq4W<2=?OKW3U@YY({0?>Mp6W1 zLO66oMg!Q#fQ+3tEu%pS2t)&o(6Pv0y!<{wKv$S9_zc7~qO%BN(Vs%0S!TtScb1(t za(~b?Un5_unLOv4$jdR$17N6t0P?7a*F4EG{pPx2!H<5gN(OL)j*@BFo7IyFH?ft( z(TfIG_HlODy(anZo8i%e(Hhe72do=(90edYt=rq%y9fpd5F0tl9tMS)f=nB+FbABh z;b_@z3c3&!z`kuAx|z)Y+ufa{ZbWVfGq@E#0K}ATNvk5ht@WrvpB`|NH^NTrkO%}- zbcqmxoJ|N3HX~61Xw}``XzxHoIt9+=NN{=Jy=RWh$4ZKzmSt5Xn z7v-zsbyD$n7#bbJ@hafDHacn4pD6e}TdTr&ZoUlZBPP!AV2$G8Oe!OecsrO5GZT%_z&jGPF z7f?iW`q7Qt)+;+g=KI{=Q}`&;ds(x@!>mW`yrT;mj^v1atSv5CE#6lQXm^(5aj^qZ z=S(SSSIuiN4J$E5-udP$OK&pxBUyvaO0TsR^U5~(=2UHYE}fQp+gP`)NofVJ3Nafgsyr`L=a-_b$O?-xcD$vC)QC*}iB7x?n#o zP6RL1<38I_Zg7vX`5WFgB}>3fKgX8TVZ0Qp_}mS*W9%0}DT8}iv%_}y_H zcPE8OQIouifo&j6NY(sd#bDhlrStK9&-aNr94RFo%2VUG7l0vKa2|3cDR%@WO_n?t z?zmZP!m|Owm;*o#Ga!qAiK=>C$3@OFbkb3j)Cg$=a)k9SMZbH#OjCNw)8W8-Is8p;+C1ObGmTwi@S-?jlON2 z;YgD>DkR4;gG2vgnCzw9=#$;^Ne-$Lf}dgV)zHe?hO@`MJvebSUu~vGg;3Lo9E^%FZl2LF1{t5qIKbSG8$+ zd64)%RoskomRqmC&X%s~j)85YP!!0_*yDNtEAPo%65%dl%>7~hm>D0EV7@U$G#nMSyz{Qjbdae*iebO3M*53j$eu)R!3 zRNOonw8`f&tK|}3y1yQ5*VLzDOIkOsR9M^`G}UkkT`e3HNB*aY@Fl-;z=A+*fEG>s z0kOKWC~=b1YUz|v2IU!CgI5m-HoH6D-oyiQTD7=~)|wEY3_hElmFf3kNmeN=6uF0S zfxPC+)OwQO-Ob>7X#727FE1w#Z;1hu$7|w8(eZqP-`C~oI(B>_3A(elUcIkv@!n4@ zRzVjG3^uJ}CH|Q$POX1BL^NCtkk{(ku=p}uF_ZukWw76DCFoy}3Mq&KrdQhd z>T>$bMG>{l{X2QaZV~a=3baR^b|w+qV*hBCpx(#F=^yZ=K;@rga<)fG9;uEQFx>EX zyct*!cg2X~0q+AL<9T@-LUz3uG#>P@mMaqXO8vly|5R(kRbM2 zRI>T=w|~N3uM>UR8-91}gYg_TE3prQh ziw^Trd4-u{$DA5>jE#6#=iZCgJzIp~JgXY>b9r2wU-hypxbIe`?hrbI#(}$obGWxv z?^=QSOCJp$vH`~SDkVpE)f>v>KxH;@Q-wEQC(=vMUOOSvd7JmgZ~=fB7lC~|HbG?s zV)^bVV|*44DsbCiy|KlTF~OIgnR)MN7d6R?mOh5I!i9{kmkH)Ge^i$$esROjFxS7l z{wSm>dU;zA+ZPhGV+S=gnWN|Po1i|(j;eBGI^#ZKa8mUTukFf_^)Yvf_&{lo0FM1$ zw(b1QyL-Vq!nsXGwyNDU*0XK`*n4U>iw3LnDe#EJ%ApiF)CX=Oqx;25C+n-yJ6D#T z=Ln+i0xr~yA20oR;$7ZMreh!nanGd_)oAyC56*$+`iB=eUP~$Nh~$!?<=*x%HP{YM zrfPCik=xdihX%TRV#VXyI-|r!Hn~8O{G#h_l~I`<_uM%VN;cxpq{brRNAO}mrI*!b zX`wL;sRC!pin?F0H*K(=?4>hq&EY0AlSjE|H3EaYo>z)&qyP#&x9tgUiM1*Q(aEg9 zd14%_jO$e4EwbLBNeV-}v5Fz6Fv*IDy}s<;)*(*r9_A}rzUeQB0hPu&f7G3JB3eJw z{Ofn>37Or>AwC|1_r2cNe;OsegeO*4*VL^INtOZ1jKIU-hgpWMf53}K70Q>zCvv+V zHn1lbJ^4!h`-c{zRg==p1vm2;(3${(N^=%8##cYn1p5I$FzF*<-T3r;AX3Oq8*Qiu z#5_uoCY$E}c)c*^0)hK{IJS}#%f}eSgC8QCLJEMN@$1yg-)-adP2Vn7$7o(>St{uE zoyC1%6&Z`Z@^%*Y*vc>b??OS{jX}DvnUU&DdhX2WcQYd|SY@Xzxp$o~d_+NXgTL>F z>Z5~mPqmXlLv_tC1d&uBjTwXw9{3!&Y+tVjLjOXlAn*C2X&`|#X<&$$MmamJ-dpsFsW##9uJ*Sxi(Jd|s#%th3>nM@6wH1YCXwOq4 zUbRg7Z4(uL9*7Bd^7X1Y!WPh00gEDm6D8qyt297C zT(U?TzY1EZ9dw(Tz%!9#D%BS}r2z%|=9HnbUyc`_WWSkq>aiAxcD8+`O7OFWE1DZq zdm#BGR@E@N`RJ_a+-2ANLj8-5b+dNoO55}IlckjPnTd~Yk1-joH!W7=E5ar~fG#ORfojXI>?EMv z@rxY}N9x9M`c(T8d|!_#RCZ%F(k1nOflEUap|Xynnyq2Ba}DPW#XKWLch~Lm3u(CS zqQZ{rpKlN7)&l3J{8$1K6_tySrx_nD>cS8o7XS!tsC@~W7yu^lXvU@fAi3Lo{x>|e zr$67Z%+VNW=M*~UYZc7<2o{-z{os$1e{WxxWBr(m^|1>@ReQql1M)>?i@Ue=G+kJZ zvRtb9iJpC<`^!&j%jg4B)Wfb^t-;`3{S(=)s=A#wWgwC)-%pnvzE9&|tG4Gzf8|NJ z7yy4fy(iB<+S$T;)>tSoyZo+fFHb$LlU^J{> zjkz-rorGs~9JlJK#7ZHS-}bET)+ky~=C#2HH=f60L(!AU_m=c#BcX|{WuweAMAf^g zO4GCR3`qaW;XbCZ6_nP+C8iR9OmxxTEJ|2j`7wjqE2)q&4G`vXnS+pGw25CCE`EWFHKSrPvQLtnf80UlIj(J{ zo){Xnc-o9LFhFusgf^sr&I6A+uw~z!1so57RE**=Jnw(|r$Wn)m1Aws=VK+G7~@xn zRu1%uMEw3K1r@lnh{gm^;mOAi%RSoCq)QSI9Nb`F2eb~UR??i~8{1g6^SLj6ulB8k zrmxne(l`6)%3PMqRR;W0;OTvBp6GBmLpOpKsDW(cv;@zFEf8en!fS^{;wC-0^ zU?aBOuW6-`h$N%NC;J$rl<(JG-l!_|B1Zn->x8{_2vlyR?Ak@1*d{0`F~-%{)s4{= zn6%OY@0BkwuX%+v6z79`BM?6Bl_N8Bb0m7j9TJ>@*77=Z|0fuLM{LjhDzo^Y*b;3= zerG!sTl@k~tL|kynR*Nme{@JQD>C~So(*)e5`cq}dfa{O2=xv`e_=aZ4Gl??zzyi= zy5=LMJzxfJ z6$@oA5@(}s(WerHHY1ey>bAFCDb|nT;ER9RN`|Gui~UKp z`%9kTsn9P=3kUMyD9lhFy|>ilOBgj(`Jr(`!@|S$IK{*kkUwx0GiSizlkfMt9=3ix ze5TDS&|N1N+90m7i3rOqPW2cJA9iIpRSy0MDw8=?po1h6B$#>6L)<^J_lflJBOO^x zOdqLUp;g-3)@~oE7~Btw+>kxb&Ak0Q9fXx-43zsvKc>wBhxu?mdVYV0IiuU~+R7a? zkMIFfk~x4HfYx3s2^~Ez?Bm^kP!vcQK8}@fGiW zKOQC(>6WcbQt7kDIh^?beLqG6EY5YZmwGSpFhs{7#4B4%`Ya>Dwzwq0Y^zZX*oVUu#U z1|WJRB}W-sfJ(YQQe@h@Be*1DJy--Q6Ey|L%TSDQ=V;8xQpAWgDdx3Yn%Pi@U!y zNUHr%&PLR#r92%nAn8y|HivQq{HDb&CL|8`YR&a~C-XL4ppl z@BW8F)jSUhl}q!46J}isbTJ1Uj}l8Vo+s-stbU4}rO-$e&TBldB?nzB4-trWfT9(vWSn)%Dk3202!9AYxtCZ3+E}4kf-E$onHJPwrDkP%U34j1t!=r! zkwafJVkiNLIHQ7baS*%In1hvyX3}x*m-%qmP`yK#-k*nI_emkNZ9Qu#qrPnyfyJfV zCtGog3_7g5uQ+;ZIg%O_D2e;D^lW*brK^_I+~}YaHP)#ZovOAQi-H>@6g!=Sb7%GU z=&TE#WLOIrm7*@kF`4RV|}Q zedylK-Q+D0zN}Lzy6VOdXW10!2}XFYT{%GZC&M1MZd*}6WzaunzD1$G!i5BJtZNw5 z&j*bc^g(bTBShaMpSauE4hP!~gWg30f+hgcN!7lZvkqqHd~eV%~8Tt6YtHv52m%uFe6^D^dM zd+E_ex1MWReZpC2yi?aL)A=bju72k^9BBtIryu>t*sdu=9}J+*voMj6z8Q`2xYA#Ql!yGQp|`}odso ztCRiGlf@qVm)Uuu!Ere!**StkHmcs#$PFa9*RE<0<_#g0Lqi7O@O1e}?dYHPjYm9) z1O$2mS?4Bnt&bs#o5GXaJyTNucj#HSXp5uQ(ULI&6g$ApQcze!Uq7~N*|L4kn}|AF zQ*?}4Qh{H2P7|j|aIc9p(7($Aeg0SMtf8yIzUkf~0;LtFskpk;C|S7{5A z&&;r+Fe|F#xnE6G`z(HT{!x)LBz^F|55ehXk+W546yWAHkjtWwxfMvES90wqNC?Uv zwUYr@1~LRSBbea#q?-}4)%rFpmcb$muRt0u+~ft;FPzu{0b3qqI-qaBxR|)#UeE1Y zs8co*^*_LOG34030009rJ_7_mc_p0IBfSbCfuAQBOS4P*C1jxI0AWu-Z(UrncxnlX zP$v)$A-hT|dHnnk*5eBzs)W&iZbuK9?P5W>l(eF1Z=GMfA%-%Js$$ZCUcvdy+9sHtj1wB~ zVj+eHF4`})VbnB2ARe!t9G{1m9bJQ->Fw_iyiFpEgA@cc1^4w51}@nPDBrXbq|hGR zquHK#5?c-W+iwJrC&#CoinTA?=@%<5teoqhN1Z+lv|WIq5T#|5;2R~6V5(xAHVk8U z;)&UrE=B?RVRW_Gh#>2mO8iBbND9c@I0{vHVVaZqzIk?(DYY!<95Ky z=|tV|6ny9*8}kVi=b?>6`b(IX2q1lHk0NyF7T+0996QX}9f_Sa;@qN%pZ+jp=3uy7 zw)-tG97by1dhN}WY3c^mY#Gq0h+GWTGwvR{F!N%|2Cgh({S$m1w6(U@1V_hoK)xX< zFOna%P;-^$d|*}LBlX|<#lg?a=_;m&wWRHS@bDb8@4)>R*Pd61!`Tn8FT7|yZk_T@ zRsZIPg_s`~DZYPYD9=x51F7tROaNPqy%6`HaFk}HLHg)*y-(Gh!|MTeB%6Qa#{I1R z_)*@%<*x(qLI?JkHn+d&DLfC0jT;#i=pArWgl8yRw0{A1*npd;$mWj-z4gEv_&rgA z8UbSU1CDi!%@5VG7`%5U5IL}8HJ9ApY^hRiUqdl22xYjbK7+2p9yHce{lckHajWTx ze}G;8q^`hl%8tRlZbKw&paRHCbSbpC1#;n6l7|!Lk@E#I%zo_=)4h!uRE*OxD4J}0 z@BP?w zCQWNhe;?~j^<@fQs?5^2IiZMpy1NB_Q0@#$LU{vNO0#J3U=Yw@LpuZc!Unk`?X^Z4 zcTHy5Drw`3#?ez&_wA-syh6JqO4qLR7IeNXena`Tkli&Ad9$EUqDJ(E^v-H~w{q49 z{l~q!c%fGpWZcADx+Px_KMZKw$WM%q1Pw1Kyl~w(C*-0*#^%BhH|8IRrT~-?+)rv9 zU#=bn?J)YlHR}e=3ipdCt(IuKCHAJt5$clWdfG=)kNNZ$i2cIi<9p@c{)w2N086N& z>k;TgCI?AG-2*K73@rlUW{js~D;nEiRS{VLx0hFbFA` zb|c>_@G0Lb;eefeJ>nWMKY?@=7ez!`vMQ};fP9XXZAvMS=tvnaw3wF(EfAPQtn_Tv z_nwP?o(#_}yLsxHNuSZQc6Z4Ku&NC zIyssOmm9%2tH1M277&}@dEC8@zyL-AX5e4q=>58pY5b<=EnFGYyDvljQV&t(zoPj> zZxspie8gSbzmUD0jBkPcPNbeZE|T`Uz)9$@vbYO9XPK8gBP#333;Rwu8<7;Vw_mo8 zA8rywwjYa3_dKBSh^_vd!OT>TESD=10TTZ&F_zv)5{9Yg57`kI~gkgxY?R@@u2w2c^3ibEWO(v)^J)XRci|K93Ue z=(iVZD7t_mO@78h)_90rZ95t@UuoSTLW%|}HpJ&C3Hsh3T*&yv7&J)>cC%^`4AJ46 z-?#`&^M{gfNKFtu>ms|j3Q6EJDg^Nx{s6w#e9t!?S@QZs35I&IoG&^C4-ykFnxzr2 zqCg>4R}(-m8BbfqU{*3>R*iU_O>i0Sh1JxwsWDIQ9HB0gDh-RS!?t9)I)h@7& zphT8i^1Ndw!%cOYGhT)>p3IBsYEya>=)@}V{XxodDy<}CNUbdhxMdcUW)S9i*V;!k z2#2^@FvAeAwPsC@^i`t8kY5AyXOtyIC~>x_x=Xz-Jc*{ghOTAM_23>Is$(#IYk#5P z#>6Ey?>)T!5%eCl!R_B6f$zVw^KYXkh2xt>&{sJB?rub#tfxwTR`sU+{Bz#{--33wYx(0WdJmddft>|;#ibB1*>M8BPfb##>FMBQ zID%ngqI(hbn56~bzzxgpo!g3%WGCO}-_ToG zZFHrb>yk6`BN};?3UYL2CYNY_kPM5IWm~CZojUKd@a_2yKZ3A`?N^YMySD2(satHS z>?Z*ZoiRv&*x({&eT8>f=F@R_M$)UKLaMrT*$o3v_FE#7jwRaeN_SP~EDJUK)sxLN z8qSxH&d&P%;ya>Fy+tU*IZT(8{qvrs^WA-Yr5%uabkdOjSf7jMxR>ljcigj!^;Eg2 z$%oe)e~{egb2pWN-p;4O-A`*3N{_SnjWAx=%#OWs)Dv(CRG|NaTj3IvnYIuRF^bWw z`KPf)bN>nT&gLAD%8Nl=W%~G>4p*wMPqd8f9cqAB){#6rx7<9IYcJn?r{NA^xL+5>;GNu_+DQoEuh(TISxX`S`i^E-J{M9vE6m5ycqpEzRse7VN=gUm-m+!V>rq<0PAri47N67a6Xlr%Or^L~CH0Ud* zz^0TVAIse}1C`AaC8%WPp%zZ;c3#Nrn=hyM0*Lb3vR@lTX?7jo@UCI@>rNGy8NvCi zB1(R9cGNGtLE<8HWGxL8HyVNq&um8;0c3*i=5ffjG<98m0)qDfMm3XMNwW7>6otMz z$7++0fY>2oaMmax86)j#tRa-e5bAd_O;pV2?>gIOHgDX-f93h*5U>46OpG6efA_tG zsxZ4^2H&xGj7f&PJQg~{9^E(xdVsAyj$0k&zPJ2&T&fM<=2#XFvhxkbzZ^Ye*h}!l z-A{Rj2**KVXydYQXOF--XfUfR>Y&v4JE`&*aw1?4aP8yB!Dsn5Tnhd4wMhzLOF+5x zE}jhiO`C0`xUL@g4%pt=8m;{`mBd@bPdpv?OOCoAnIn9y8kRPZIM&)iZ43gSXL*Hp zyTC09QBMkd{0yEIi8qsWpd${PK3-Q4cW4(l;k#%#6?Z*RW;J)L5PGhTX_u_k>}~kP z(?NzCE0bdIUKfCL;P!A^Wc-Ugm#S{<6TSP!BG>_Vlv-b$^2d*7i4~cBS~((*^EDf( zq4`OYS(5*uLW5j!3)`MVH|MvTGm+@3uY2j%84b@Wm2bTLQ~vn+jT(gw>~l4_6!`4* zdn(h*SAA}53vnLbf0Epr2t*R@ewP+CE&TxygA3O?1 z>r8UG(|ET=9-|+3eP(QFs1hjgTG=kj%{4bU1ETOXVd|@{nHY zZ`C~_V?#i(b-x7&KJ+)TE9 zz?2+)=p|P=Q0aNzWO%9y^Wbw&oJ1)*;mu3B;~&uEYJ!s;!V6bYk8%I}syz>d06k!v zNfjmi2t@qhlN8;8Exvd$u;FH^T( z4iacUnPfAX{U(pO;4fuHuAGs&nt66K!FR{>kF3~F-sEr)%$7rN0Jd5a-B|DjM3D2W z3Wp{60@Z@i(P&o?L(l}TjyXTcd!A|5Vd;3_DCv!vrC5^+C?c|n>9H=NqlSDQO2D!ESRZsn=J+3LOlUJ978%m>jPbkci;Fpm)Ttwy{i6EPn-KUx@%0Mp}G9Ja-JtQ$b?o z+BdG!ajIX z=J%X$3RTP0dWNTd`DpXs{Y@b6+3KQdpg)e?Aco*f4dxX$DLxo%K&l2^ltK@P32|xA zx@KB7ux;d&qm+*Cuh)_gFy_TTJ8sL^*?}gCsr&+_&ce?0*_*N$7xkue=9e4U8Yf)~ zdW+SdnX+zMD~-328uPo{=qxzAw5hw0*BZ?DvQmEe;5L0}{?i+8DaQC-F#=A=Xwf6G z&$KK=0u7y?63L(!QD(=gZDy6L!1mkg7k_DPOuHQ}Ex$)vPs4rg2Q1%Dd^aulev*zl zIj-Vq+u!1OOgdOA8H)6G-=X>dE1gQ~L0o2{ux5 z;Sy4f18)2Yd^V9&gUswQ{a}MaoYSFy@w(-YtfZ-m1n}DohQ6l+b7rj|6v?OT%M}#2 zYzO93dg4Or;-emx=sm)0>xPJ6XrusZz{hq%CG4oUs+zuqBbi*@>Wk~PbrH#P_U zyt>i0bkmGRUVkQ0=!ebBX~QY|6{8*WSM~I8e^o`lSL_1XSzgm#9M?0NsUv_>Yt$(eG%*tp$<3bi1xNC{d zar?@-^Eqo9n6Oup)8m|}$d6x<+&fqDVUi$L=l=b0a8LZ*Fz@s5T8-cn54X-qDmL8^ zOkCrRS+DHViO7?v;I)JVP}X51_)- zJoI!JAuCreR?Me-L&zfOc~v0_lN5xNZkmyn*Kg0P#Rx{Wp@{K27+mS;UL;7lIXLAOB$_qgOmj<7CsLDy_oy`nJ#F~&7&~B;kD`BDVV6G)1ABm&REHA- z{`RUy$rx~&6hN96tZsBnA+YuC`V@jmfF-oCQD=>QZenYEIF$LPz~eRa4>guI+AwhI ze$+Q%{7d~w**}%|@J(j|UVL8>1_=?s##i2?v)gaE77y-PjOLYZRnR=Gv1DBNu+a9i zYBM=$tJH|0B{%ht!c=C*^}jy~?F~E6*~1jy<@JTr{zRa>GohV2FMnb54omn1NZ%ko zXo2!E3}^%dGeV*Y08WHi?HY}{Exbr;#vF5AGUQkGza7AaN9z(NJE=c^1~p~fw<^?k zJv%8?SC89LZSK)fR8s;;J9Uo5Ly?<`E0n!)>~L0)?|6y$4d334x*C&uMJCVw-rY$H zhV9ppZPW@&UysL@-aH%0{<*XMLTI_N53IHj{(ER>C}}D_4loN^!7g&fm0DWNIOch(AQiaRS|uHE2wen% z>`^Bo-kwKpzxg-EdUoe;)E#eZx=qDX$9Oy9K8-#i#Te&j@=7%8e{Mv%eOA~>pnZw_ zVemEs7xEU63G>+w@gswgbd6l%YdEa&PXEq5CKjEGL&FLoWh-C20+ogVcwbLm^~6Jf z&nCYSNev^_ObHpVh)Av}z!*T&s~r=|FKG7pf7(0Cu&BEC-;W?6NGT;9f*>Fw9nzwR zG}16AqNGT72m&e+O6kzuAX0*q(jq0LAR!&nF*N7i{Ql=$r{0|R=ec;ih?zZm@3r5Ktb)H}5?5yhbg-#24QEZ%`=iWsT|B?O6UqKwZr2gG?kC zBrp!-{hiXg6SzpLDG!2dgu7AAM`%sF{h{h_F#UBH%)_-1nI>pD(JV*Kqpz(Ewdc&K z7qt5AX-48-K@(u8vE)ruV@EKSs_+!U{_lImc?{~-lv2``r&C_d(s8*mf#M@ZbYJFM z^$-9VuA4QWPD6}WTRa&_U))vKjBf#T#IFjIeq|L95CLHm%uNMr*ghe*5CtLb7$;MV ztZ>ZTnIGKm?s0hI{zIl<5(Mi)_*F^TsTLI05?J{2Va{+HA`yQ`Qv%ZB@`Z(;vTgjT zBrh!Kg}H!h1`&P^0rfLH(h$4jMj)fYA|el=oxpo3k5ARFT~cAj&>|dg5MJE!{!851 zs4l4cw5a~rQwn6x)}M;OJwXQ39r_au+%E{B%(PMg6L>eF@^cyRMj%Lsz3Nx+_BCU6kT?U{Y0wrX9ReLj(mxM@H@s*KG#={2 zB8nTT+WmEHhe~d^{w_;{GbV-#fy}42m#zjgGbFYN9oO(p0c|mQO5z#i`fx9;q9PX1 zvicE=I$`7kgKesv2pClMw9La2$uol5H!4ag30gr5w^u5)MWDmTYlhjO@bu1Qgh8QZ ze2oY9zfNo$xU5Yt(_-cu2u`%3st4pk7&AFjVS__T8i2&#fB#JkkkEkN0m-*^?~tCH z_VB4^#Mx}~&S4syQ9$y%7z&BFmKOeJ2;K(y%<)Vf4H?!IdVo;Iv{u*w>9B!EtVf3! z2s%Ip9IbI1t|uTE8}AhiRVUv8z%&JL+mifv6{2X)n97;lNBXUx zQc6p?Z92WB1?A(A2fWv3t7!zsE^vd){hZ%L6x2JyY;u+ueJIN~xDC1GI;c_#Ma7id zRWeaNq>iyLU=t$_l0Q7Ue-?CYI)$|!jJI09KQ8r%P+9Pu%c@h!q)b52NQVCjozsO? z?AFT)sZO<*hfX9}kWFD$?_tq(M+&gKYVyME1^f>C5Wc{y0|P;vx+bG`2tWstQO2&nEg)az86zKGxKN$c*Zd;ZdsG}MBO82}L$kC4 zIj-wkgQ`dOyZgr%D=hMorkOZ6=S_^+n&?D*Ifxtg|7__}9Ht za^z|X)!l34u3H#FMq(7Q-X(a2Jg}nzt6viBhYK z-b6Q|v&z$1w5?soeT7N?jb(A-qNJHb@K3haI(n>aP^i<&`-QtMrIV0_i4oZQ^5U2v zSe{isHVE_*mO`%~aanEOvjmk5Obr?h^5cw&d-V(M<)c8Wl*-}O^+?sCDu;5oOvo5O ztW*F?RU`ivyY(0=M*uAf@1R7W}&txC{1OLaM#lFOTsT zzJQZL-|L$BG;&>q!MVx7`UYy^{hBC>AB6z{!O!l3)^`U?TB1DXHX7+$JYO?>2LJ2M z<#(I1;EaGELbKyuvVliDP=*_?OB{I&)-0GY6ps)R$Rz@RB%ce?VQY1t|07Dj@!`k8 zj3U~W@?M+VYgnH8a``zI`(%kobaW^`w_PObi~HoQs-tfkYGE*0E#c(T3@b3gC%)(j zVFOxxx3(g)eo8g$@VKODXX?I?E`l9|>vbf$l=qP8QFWu-Kh)3D*db?RKnlDBzMz0oG#fv(0_!DjXP`S4 z9(x0L(j*i$0U-%!Pc|DU^xuZq7`OHN{=>G`EhDmA52LzfGU0OwY4P(XBxZkHfKC4< zk5nI1c`moEI$gHu#ww~GJZt9TvkoGazT)N|VB1H{TjeLtsczY(yh;S?72H$s6*~#; zD%Ony^@@^+l{ysu-SdsOJ~+T&FF{?E8FZ+9O(WZvxKb zJ%^9)zkDh93;+RM3}3ubqe*OM{}4SKiz#ZZdDY9sO!>Zt#_lx!{cS)G4Mp6J*5%Wi z9ftu_G0QHC8CRyES_gwwl!A{1y5Ja%)J-Z-Tl1~Zm^^{NImF08BC1ATB!Wfh%cl*^c z%(~m2of5pJV$2|yQmnpO`V4+9%A!bu=mIYUFBm>G@J+E+E!=Ta0M^?yS4C1B*wKYn7jE!f1;B0`htH)vb@_{9 zJ#~&wp9gy`M@A+|^@{84B}r+r21FLRY5+S%Ys{3}XF zX5Li%@u#j_R(ZFh`hL@I|Izwtx}!3QUZIspjk8DU%z?*uOF#E`%st%5;Gv)$ma>$$ zA^M%g2Yv`_H{Sg1opwn1gn<3785JWG2pveH9v1N20w(#1Ug<{fg6th%PX(2_S^h*h zJmFHt3CC|E-98^`wcYzr7!{lR5)w#qHJT&CYRJo<@svl(yQvyKNp0xBB0m?GvV4HMWZ1Iudb&V>)1=K}i> zihT_%`)ZDv*WGJwu?oQiW{RK@PJYuaWWQ_j4#d+ZC3ETc+|gNHowe`AEcJcoOVM1T zOERoH&VE;&H5$0mpI-vn>;u}mAZB-LR(Ymi|I)*P%HL6HlZ{lBWF+>6RdIr+aY%FFdTMp>s zAEcOvV@*34?8n|doQaAFDG(j>V?CiX8sS9$H6FV%8>VyXBnJT#i5p5*)Ts3d0bZTy zdIRybv@t21y_{sXG6ywzVTmlECNWKH_V|LSm^6EjW*E=SXBrW5Ik8vAf6}!&Y>4_XC_YwYML8Y8&tVQT zuUn@)h^L#7IyMIivAtBxCHaFM8Jki1Gzq9^le;*b4TR>q47uH8iBZ=2T#VoE=KM^T&4eB;QDb*poMmzv6cVHw zIRrt#5y-Yn5;IDn)<}lQ4U!Bv}e>} zfOP6c4umip-VA%nXcG3h@KwCW1eeS*|ghwPl(L?yDf z&h{Hikf(z^MCU4=Hpcv(Nh_`Sr?de=aDHIdPLp20#iMG(<=;M}sI1hO__G9zY4_Jn z`c)}emn(#-jcjscSw!=Bbi69D3V!JRCXNsHzI`fqx)67pOoOOO{Y^~NApWhK%jmuD zdK%i=fr$T##Fda@K!ei94po28hLYaZMc##!M*iB*it59zn(YiU< zkF45zTvn9OC&xGB=W8Zqh*!wYq6nmi3pri# z14(bntAbVfPGd`!&~%Nhay7;Tx}#aYgLa>M#P~ui0d$1r*;WkEQy`H~CB7CP<39#_ z+XiglOF0j*1kgi5w<*ivMpbF$B6KDB{$;n_jrT6kma?eJ1Z?4WcxvFT=&LWimZN~_ zA&jIouSUJoLMO=MMptN{(sv}^-44Ba7EkAK5=9BTqC{&0R>&`Yph>fGVL-5uH1p0- zUcuYXEy+L%xU90@_SqU3(NI_uC=;7=$0N~wU%V~dRm;+;Uu{c+-z9$KJITz=aTFUU zR09canY}3~8d{JA9tKNP`^&c}kfHyA6D{9Ug9~=x`471}m#!-EQZM!R7GQz23pDo$ zOp}(p$DC|Z$(WtG#49qrKLM~<|1LVcfO2(E?g6X(3ns5TZdR4+!=1|$cS;o{!vn>g zW|5eWeF;R7>f6PY+}`Ln!}YlI9}Eipa}zpAx|nFq+Q`%zU0Fr%unTiKQ%gi~ed$ARF`=>7JxRta|Vq!Fa+s(|cJp zsq!jp9jE+hpz5At^qWPX=57koXW%oy4!>k8vLyMYW=gD(bUd-xAR9gQH;gmxRNw`q ze<+Eb6IRtu;DI`8sh`G%u@+Q5(YIsKxkH3i2iS(Ac@q#Zysho#ZTvbHrWbW-iE*E- zNe-7(Px}%NHv{2bN)gfd9$%GBrh&bWGTOOZ+?gOq6v& zLliLrNiGP0eAg#CdWfI+>-D@New0W-6cl-nGu6Y15G=>ja{MosX(y|KJw;81LN}PS zpHaM;?%D-pEeN{J40L;Cl9xLrw{L09Vc{^r4;m_#)8^iSM{8bSOvyo>W*7~n*l9fj z2F?YtEGYLNIgCzud6m9I`un{nncbgwVK_l@x4+dUhl57~(VuR{0o}!FQbay0f`F03 zr==B5VDRS-i)Q1jKjyd!#&^!1y$$uKi>mj2z+4y5wR-*YHW%%AwhOOzUC*!bLI4)< z40kRoo{9Q7x$fR24#xfb)-}OpM}vFm0g#f8$a6w>p` z75;uu!1D|OmoLE@KN~!l(_S}RYG?n>w3HMK)c!|^8W?`P3jlmjyY~bSI=A0kSzcwt zQaqeV=@XWx)B(9O>^V=Zo+`9vzb%HQ0thVZ$(c4Q2^wf$lA^GnZ1r*y%HEU?v<_Z< z9P0%5w=mC}lE+{(zQ5r^y-y;}*r+s^?pyO5&cLPgodDnj5VpVr_9ghX;9VEbC6P}| zpt5gT)z0ev1}p%_iuNb%K4cGmwRtvp6@xJsS*;vm5nnO=M`mvu^y*aG#vr!(%I47N zWN*};;5x#Vk!(O_XQ{ervrtzNQZxSUAg*yCt-JTvdZ)$voj>XxBjKjZ=KlnNK|?59 zQo@_Bf=>_6HdNce7YIC{dM2$UWZztagcMJ~vF~!;0ti1Kc$L;Cj^_yEP|Q$OLVU~p za7DJS!hUaN1O|z{>EaA_J8lrCY9`;b$PV~qOxvONOt&%$h+r$;BYT-&vsvzxpY=0% z%v$@{L{0Xu8-x2nj)!Jty7M+>+QqSY{MedK6N?6zz~k=sm{QZO^g~auw|hSq?4A61 zm)RlEN=!O^wDGok0eb*HI&}%#KB{e7`u9QmR#>RF^d-Cetd40H-w~>JA|5|u5am2- z!geg6{UQe4og_|$(HTiwiQ^|Ry`DEaYCqDI(LQ~kypd=$7Rz*RdL&k0utH4~T4_*w zRCu|(NB>5UJdiV-%{bG5bkT}`TCJ=~-8qqyE*2Ke!7K{l5A3Yl+dkAw|BzrBUE;W^0B`y)iJ!zJ7)HVA1Zo?{;=yG1#CMZNlhp&dy_Jn8bz=^Olew;d%Z<6BOSBFykcy{$z?SnMT1orEX`>WuE z&pEI+Ux48LnTY*!U8BJEi}Z@*mjHHw1X@bXD8PW9Xt|_a@@T>FtaXSPjJi`7@v3yj zP(kA#yijtuW$^yX(idNhWp53;JgkQs2y*=fxWHzYAf*UIxBs}?^4ei{t=-#Sr_ae$ zY1lJNh1vG-Hm%^I@jW=1;U2nsNj?5>5|jIjt@MOpbZ4GyQorbu)SHgSIruFoLoAgI z{!tl?#<&>5SFkKHNb@?0vxwcj*C$XuOt=0aDuNpgo@_Ps>e0uy+moetu%(4w3L`k= ztL@Y0NWjMZ3J)3k^)jWy(QP5X^M04pw0_ag>TX0Yv5NqIC&H6T~F<6q|lAEwDF!{8uS=#~d#Iv%@U*)&oBn8^Cl( z!$E%i`C}#Lbdk#)&$-;}9BS}qh@y0^tADtlSBuW0TFeET1%>OJA_T@dce17IU>@04 zP9M*_ms+?BlZqmgriVgY-vKK_tgWE9us&h!NRWi1(3@_@62s-%Zqewy zn>MO8*xu>STG}cawzInDa?86!8yca&1_Ys3hib8Ll9k1RNC_pYVvIQRf%V5 zw#1|r=nvW}`Imfh>}j_xmh2$2J)}pJ{9w;FPTqGo|I5BqY@R@R7hcVTSqaUX$b$ZM zqo=-ug!XSkgJoqyE_Zm_0qJyQ6Zr+Df^A0Y*yq@{lZ!B+# z=P3LAt*;FhHNuP&|N1K73Q=&`@E=ayXi;(Xu30|VoEqN$G9HmtZAZOaN|svq?n
  • krS+*6$;OY6g|qYas&br7UUWzC-B)A_Uz|_BOuV zBcTD7?yYV$Ue~hAhhTv%@m%(#a6hs55ClW_>iAaS+wLoVJH+C%i{(r3q;4d1Yv_G; zXweaikGDCPcWR}zF#-$LqXFk35dRI+iV%G@5pTF?nX8Skwj$F$t?7eLqSWsUoO>SGx3C1DbYju!}k7a2C86S5LrNw1<|RN(<=Jf zoZRMy5VD-g^clzFW(rE!#hEUdPrsDP@Fg5b%DbuKSLB+PJ-`8DsRd=H@(| z%9-tDo-2TRptD$v(Es-p{awz?lH~Z)kAPX>u=vUCf!-Sb+nf7~`YQ1U{yoweqf*ye_KEO3=2#g z4%3d2Z#GwzR6FMrV83cyw{A^I2;Eq5AL3pe|NKMVb6LXS2qrAW$5;_;P}*aojoLM@ z!l$_WYpK|)%iF3SZur22;97JOcCC5-+Z27ld@|;t7m=+(`t9ldr?68|9$8Fp&h`^} zViofXO6_4%PgP=P z@&Ke`h6fkEUU;N-ZmV*^l4<&gYhtm@LBHGD=Hm+}e*|+CA3LuKpe6U~8bpv(5mH@x zC*D&1TB2KDKt^yzNMAlv*(;Du!;fLm0d7jD9NzSE6u3oSFX)Ry;wrL$5+@{lR&ESP zc%;lBCEg77XPXA8b)}&H3y5ju%R8^uW1iFt4IZpTZB4=oU6(@=&|i`)&~jH*6TGG+ zt5E~?{@&r%ooaCD*P)~@r~O+btHLv=p~9)HQOfzUQYqmL1QIaSxIoH`w060Q2mJ~h zp}*GB{c0X^0BvPiRd^NUR_rJ4;avx1mS&V85+~>h#bvFwU#=|*|6~6@$&8^%Af#G! zNna0x{o*+7-YC`A78owN$dajzkvLx_^^mR5rf>oZjz0ai%@0lPO}|}PHO-Eup9w39 zpx$%)AYt~hhGcn&X>TP8U0oga$M}Ai+mjVHX;1~v4O^n9z*%^+eGLg-r&0T>_Pwye zJe;?>iqe*-NsYGBx6b2BF~r#u0|1O|dK%{KU~HFG;KfC^X(({Y_jtq)Rqy3LT_nN- zA^(g+Om`#zVTc0&PCYP){gsS!u}{e8tY*`Cz#TSw9_%YnXiXmD1v8M9Xey%HwvZXpom1;INF4Qc9+&t zz8Kor(Ml&(u;=RB?f=%ZtwHi0F4x_=HDC2zAmzgmS2aG;+CdjaP#%CczG3_5W^h;mP2*X&+pJGl>^c+wE!5a0YBOzA zWg{D7k?7VorTy-reA~$WvI6^WDup^Tc64>JZC|YJCVMde4-1zoVL2aP;-{vIHBMr` z+XW5*(_*X{H`wA4e!7>;tB{s99Dw%Z`}O^6WOBtI7OJkBx8H{2BXVA7?(jqTGFpZa zps7fd?Kx5O;q;er#!j?ex*0EejOPp$W!NjdhZ`^O!DzD)APYD2kM+ZP6&~w!g2f58 zOUQeFN`JBJnqqXS!Og6FFX=+NIX4nJ1DN@=K3%JxGea8{h*7a{S!loXz+Yq|rAXVa zo5prGGymY7zbOBO4deJ4Q@$;2iFP;Am}--{@d4+0ArCF>J>R~wHHxFfe0vq$6IGY>NfF4G6$SyD39WAMq|B*A^=%6GEuO0`YFyv%-5vNPX39*rT|JM1v>R0W zFQxZ#or0ay78o)x+48@&(7eMDU95Ar>M9)zBRsnY4ofHL5MwgQnd5j;^(h9dl?CBs zRq|mlhTnUQupuJ>)^kacX=A>@ch8IIKd0Yqcs}IwAw*bCL{U#D zgz&@$AGG^)MrCVN8-6TcwtpG46J=97Z#d1G(l?e7Q&x|0lP{Mtun*`2iq_{9?lHxz zXt!bAV-r3WZ9AnScU;&E)w{SmJj<`cmQLe=fs-g%{xgEP@~13T`UMA;zq`ubcVT<% za+WXleHUJFY*0m}7m57|EcZwj>Zgs&mnUVN;6A?cN@YSTC9`^H>FqIT3@6}qwx)ac z+)O!8+Y7>0rF3Uc)ynA5E6R<)togo|pk=z4*_~wQrMSwztlj~>g`>2>qV$H1=lFu5 z;Q`WZ0jS=(y^tt}6b(Es>c@=sH&pXV5dZO`YVnNSv~!*^7oLt0Zmt1+%~JZTKcH3$ zAcLOv0VMJlMj-km7S4_-)Kn!KjoLXY*920sSr=#NYQM$W=uPfy8%(T|UBOb1?T$_x z37AsgVWl~Izyk}G@O+^h*m4D63{x%w%gdOOZstf-dW6u)Igg-IWBT;vxCJ{r5 z0(S`vSGuPOdp`OqPQ>Y@kXM+pi}QV4`~J&gn8*%Ke}jFR7Ax;d!e0zdl5J7RC1D%{ zeFw$_+vjs(becU03UCzghd#Rh^+voENI}G-i#z!yy1;5(vwnHFLA0t5JJ@tT`B$Ed z5b#-m?S>eCFal!nl}Ml;F`@KTnC0QDHA`>L1?Wv(&gQQw-t)=umWbaRvTXQb55_mx zYu^*(yK2Ip7V; zn>}dFn@iVD?my^Zq7>^b zQLZze>oBT=g+TL+Ab6Ny8xe+mx4WxowNSpK$nwr~^S8|fItZlf%zWd`FnguJ8XZ8& zryhy2H>t1+ohP{AA&=Yh{c>)!+8G9i8k z(Beakwxm2K_{Lqcvl#Zgg*#s5dO+w45GYw|zgG(OoNV;i++f?Pp3%sK@A1a$lVC$~ zV0Ysy*&w6rA1%5kcn-<-w>=$u4Y8i&%<_#7nwdid9UgvQM6Af!w4MVmH=r7<@$O3R z#Go?CN5?~EkZZirip3NTTPc`{5x-dInxz1N&$JNofN_+{ZtSir+FcGK}~E)95U%x!|#@#9?5G#z98x4&V6S>N8fXU>oR6|(z^p8r*z zjHd{5>N%I<&I7C+MsE6q1>Y!9x1e^fSqC5=fz zwF8T(5au*pF!GWuzo4BVt+8>evf_28UkaJ4kTIuw4v`QiYtxoD5_&apH^grWdA|7N z221j~bjSQ3$BRnUCmZ{bLj`X(c&(e1!pMb<_Eyf-r1%$kKxj>vz`b0UZzDE}Y@x0L zD_b;^TB9YMGbo@Tc@_p}a#K%ctZn41RwqwQ)(5NNa~w5O+Q&!I`AEnH3ER2Wuv#zw zClM)EJ9Jz5X>2R)65apk()hoyMt!oXF#-+wXSRsj)FhXMkfsvkBVuFEmYO$Kj>#OR zFAkYMtAGB>qupe2k!fyaJWzdn;kRL9SDR-b11vV)#2??3tT<_GlQRfSJQMy<^6neb zqw4iv_lvDR6$?KBAmSYS&pn~Kb2qAbHl8IehG}u?E0Q88kw==e@D9v?bxvA%*4L~QXUjk%-s7F553bWEE2UB z*&nz(`n(*@wSAyLba#Qt#l^IGe>dnmgkD0JeF``Q2W^m2f#4~z|HLsljEsIluVE*Q z8>HTXyr8WQ4o7*?3c1Nk<41zUZm$b-KkWL9op9qF!5f1>6#MYImlpCKEbi z%Q|iMgZ@30O0ILIoatHTme+1frlR#H-BwP5nIS!+Gtl`c$+t0~_jorMdkKr`H61RSl?50TDUb=Vl;c&Svm zdvB@%#0p$iB+S=xH%gimBITzR%qr`fZBB(uOg#Lw1tCfph*pMNvr2e5;Ay*kvxP7( z91HY_Z$a7o^RRG9K-T79yR>D-EPKnimQ`n!ThpfwmK7HEI7<$br<0o;vfQJhe=Zf! z?2^DooFrh09sH(;0bz=7>f@8{so(Gg!Qbw#rAIWMO>wqgknZtGvxCJGFWX9`s|tdO zpJyfe1{xu*&$@AINOyH$<7lo1Jv|!Sh7Cqyx=+5ct=DNE_aDt{>7mD87R1eo6Jy3D zw~d`>goeJi&0%i%#QsVu^E|)G)LUiwb~W|zce;6-uzcYtnE_aP>2(sC-cRT4<&GJ8 zT(^9GdMpUSwX~f(p5Ih$%~;ZnxILuan%>P(%YUPq=Ea!k-dJAw1IgoZ!7IgOGA8UM zWacwSXWH*Cv?r#M98(^QxGD?bW}qm?SkCqf73S{~zLmFbI^J!S$y2^}<3AaRc;ZN6Y9!it8*C5=3gCegj|3K>k|z=jVX!}uKkV51khC;)Cw&f4ifDjs=k@q2 z2=_wt^35wf*n&Kf_;i8_6-0i)plbRDAz#ns{XZ`0w^JYp4$Vg}5XVio6X$M=@2SQ$y?yn0#L^Ws`P@GdYONNxwxB2uQ3N zMhb$34G7PA{Lak`tYiDn83?=% zeY7r59@vV6yAay_(H3?jRw3&bViDdZuoEMuUi^t`MmRkf&momo%JmM)G>S z+5qH-nw~In^jE$r35OM=qM4|}9u7z0rq*ZzO7Auqo&6jT5#1osM$a5IkcF4{bT4|zU@T^0o zDapaAnk@)%loC$4<0X7WyJq;JbrA$Mz138$A72&~|9*R{#oo(PKFM^%m>%fYGGNwR z_mE$y?FiY(7>lde=j8*aZ9;OF>g-qrA-H_jz0+^H-j3GX|HDJ1^T_!hx9PGTm;m}< zAUWLIyJ@LXVb$>@IVe@w&XmS`<(BK_V8UkpUQF-NmRG%D)M`@h{3bU6UG%tjGJo~# z*6?gSJ*tlE?1nOw3;uDExK=y2;nClH2G2!h1s&ew(K#2tG>})}l6WIe8J;(4s6=rr ziYjqzg57Ks7`CiY?b9ZP@k>fY2JKhCzu$7;)@9iO0UgPZsh1sSp8+dnP>o9<&nSxA zu5+JkVfcD#%pSKNb(Tx7;au@dQ#-Jzoi(?9uEJM9b&=}bLc5C*7eAxNHwbU5^*_=~ z!OIo7Ranl2FaDD7PPu~Pai7f~Yi&i9NB~j&RZqWMSJ`E>LMuMiR<|V0nOqUMG3oIG z?riP4xIb8pP?l_Ru$?@*$0@9nX1oRsBdKT&*{$6s7;i_d_zf?i0P|U zqX%jgJ=^_xyL&~os4JvWN+w~Ts$&6=`*4}nV%}4Q&H>;yPJH=KD>TDVl30sW=wbqW zUc?tbdIzvLI$5elKuxYU^ZQajr^>c3z}M3m<+HB~ZHJb*&@BPG4Pu&*q9qg0WBLE1 z*_unY6_D`3&~Djw#Q!oB0K3w zQp6(u*>L-n*M0RTe{a(||B_8`-FW4Jej;lBwfQe`4>TG&UnGrlrJwX`(oL@kre_Xa zYu%ydjKZJkc#8FDCN@0A?(W&xS~UGoH)4V6h+TvBmb0_mG9Mb8Hh?7Zgv z&mbQeD^o|5o#=hx@%V;{)XC3$G3#f-S!#Q>6A_L0Os`G8jGul+GLx9UnVDtR%@=x3 zf2*${Z#vnX06&oppKX6ryMqv>7)};~+;++IA zj(7UHadN($)8eDZB*0a0aB0lcFI^EwUYR8tzqYDf)KPT>a&O>~lkzNQKXcd(Nbr6- zx+Btcpv9;rt=!nd+By|MCjHLIS00<=cj7yF6y!`cK6=Dt9WFx(cfY-TGo=wfZFYAG z?pdSv>W4#L4O22^A;J!3bPzI*$bnkWzRF$OeASs%?i{NsLt8Txb7sHO)hG7g4t^!T zw*E$ovCq}!N>4tkS^F-tcBoQ3dH5FnZP3ns>XXBR{`2B$e@@~~*71E^`nuPu(!-m= z-Hw9<)x^)iRkkJ6lW6%g{RM)mOc(RsEo5ljpUiy`TX~V z@!Mv%0(5X#+X?skoBQM>9K+iBx2Zt`c^|FA&a+?YGNE&(-gwB}fcZCbSlE_Fdg(=; zFF}6$o0#aU(DX?7ZHc2xui$3y=P8Q<347|r2x)B-zS&<;m%`S$&dD4@z z2|*(I<=a+}1~{^;MTdphjTDl$GY|!NW7F$|H*y@-FvGtq`0$xiY~F6h-sAMyEV80= zHsKdrWdsHzl$OL1SbsIUKIenLu7I0y>q}CX`~Th~*d@XVyZSYI_1W?zC@3S%V+shw zXDH$n)YbZI#gNbS+g4=G4Fnt*PG|e~Ppsu{Eky1QL*09CRkD&rAA+6KXNVZ=tMi4* zabD_5nxj+SFI5*lz0M2*ii`_^XGv!Q}MjfEI z+tL?Ea-f$2wBvvt+0C93Om;)VwN}eXH=5;qklqre67|=0NsU`P-~KO%U+5o&^ovqr zqY#a`i}Y6kkT8H9p;O9$i%CXDv9PKWV1xn4=&s5Y77n5WV`z9-*9;xk0m(QpO+5rD zEClRZJ?5e@)-U7sU2^$$eJc>>mnGN{JY()&hyu~xUm-ua0LiUb)8(eR=ruL# z(e3%O%ZK?h(#vXXYFtb2v{)wIAq>jsg+64VN#u+IiQ1NO)nB9eUKA}Aiuf#XRh^6f zv5Gk?y7hEN$wRL)ozu5DZK=8~dAus;T(Z<(>pH9w6IZ}W;CRK2j!)G!%1lrf)AIe| z|A=#%c^=umL|GvFWK5 z;?_R3v?O1G5p{79X0VzTQSAd5!<^*4(G!jOB)|aeg$F1-4@z+J9PR0mqU6PsgAD@C zf-RyfxyBi9u~f)s-x(30QDTMSJAr8h^94{+fbfAf1^oj9$e@y^FPyfX6D9fsvchKR z$fB?7`?P$Hj?FJT(8f`LwQyy)GsmutCn-Lux{F zw?KrM%R*>z-I@WZBvOV+Yi_D~*ND9$1{_f>5sw&g)B2=`5%;n;mZCuR^G1MkO7^8d zu`-n~ySWK;zj^neF9aFCm@qJ(k-xX-H-meO{0ffe0Im)E3ckP3fZVvNLe_ug_Q_z# zbuWzl#!5%fK`q>y7!<;Wk_Rgz6F9Z->Vi3NRXtS z(*J4%32b7q|LM2?)m9PweZ&Dj{osO0#H&C3{{O%9|LG2#Fnes>+T~6m_dvma%D2=M Ji{#A${u}PSD#QQ) literal 0 HcmV?d00001 diff --git a/examples/PhasorDynamics/Medium/NewEngland/CMakeLists.txt b/examples/PhasorDynamics/Medium/NewEngland/CMakeLists.txt index 5c22b422b..bd9b0fb2c 100644 --- a/examples/PhasorDynamics/Medium/NewEngland/CMakeLists.txt +++ b/examples/PhasorDynamics/Medium/NewEngland/CMakeLists.txt @@ -10,4 +10,3 @@ add_test( NAME newengland_ca COMMAND ContingencyAnalysis newengland.solver.json WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) -set_tests_properties(newengland_ca PROPERTIES WILL_FAIL TRUE) diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 515c93954..a61804613 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -102,6 +102,14 @@ target_link_libraries( GridKit::phasor_dynamics_components_dependency_tracking GridKit::testing) +add_executable(test_phasor_converter_reecb runConverterReecbTests.cpp) +target_link_libraries( + test_phasor_converter_reecb + GridKit::definitions + GridKit::phasor_dynamics_components + GridKit::phasor_dynamics_components_dependency_tracking + GridKit::testing) + add_executable(test_phasor_stabilizer_ieeest runStabilizerIeeestTests.cpp) target_link_libraries( test_phasor_stabilizer_ieeest @@ -147,6 +155,7 @@ add_test(NAME PhasorDynamicsExciterIeeet1Test COMMAND test_phasor_exciter_ieeet1 add_test(NAME PhasorDynamicsGensalTest COMMAND test_phasor_gensal) add_test(NAME PhasorDynamicsExciterSexsPtiTest COMMAND test_phasor_exciter_sexspti) add_test(NAME PhasorDynamicsConverterRegcaTest COMMAND test_phasor_converter_regca) +add_test(NAME PhasorDynamicsConverterReecbTest COMMAND test_phasor_converter_reecb) add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) add_test(NAME PhasorDynamicsGenClassicalTest COMMAND test_phasor_gen_classical) add_test(NAME PhasorDynamicsLoadTest COMMAND test_phasor_load) @@ -169,6 +178,7 @@ install( test_phasor_gensal test_phasor_exciter_sexspti test_phasor_converter_regca + test_phasor_converter_reecb test_phasor_stabilizer_ieeest test_phasor_gen_classical test_phasor_system diff --git a/tests/UnitTests/PhasorDynamics/ConverterReecbTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterReecbTests.hpp new file mode 100644 index 000000000..106037268 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ConverterReecbTests.hpp @@ -0,0 +1,535 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class ConverterReecbTests + { + public: + using RealT = typename PhasorDynamics::Component::RealT; + + static constexpr ScalarT kTol = static_cast(1.0e-8); + + TestOutcome constructionAndValidation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + PhasorDynamics::Converter::Reecb reecb(&bus, makeReecbData()); + success *= (reecb.size() + == static_cast(PhasorDynamics::Converter::ReecbInternalVariables::MAXIMUM)); + success *= (reecb.getMonitor() != nullptr); + success *= (reecb.verify() == 0); + + auto wide_band_reecb = makeReecbData(); + wide_band_reecb.parameters[PhasorDynamics::Converter::ReecbParameters::Vdip] = static_cast(-99.0); + wide_band_reecb.parameters[PhasorDynamics::Converter::ReecbParameters::Vup] = static_cast(99.0); + PhasorDynamics::Converter::Reecb wide_band_reecb_model(&bus, wide_band_reecb); + success *= (wide_band_reecb_model.verify() == 0); + + auto bad_reecb_band = makeReecbData(); + bad_reecb_band.parameters[PhasorDynamics::Converter::ReecbParameters::Vdip] = static_cast(1.2); + bad_reecb_band.parameters[PhasorDynamics::Converter::ReecbParameters::Vup] = static_cast(1.2); + PhasorDynamics::Converter::Reecb bad_reecb_band_model(&bus, bad_reecb_band); + success *= (bad_reecb_band_model.verify() > 0); + + auto bad_reecb = makeReecbData(); + bad_reecb.parameters[PhasorDynamics::Converter::ReecbParameters::Imax] = static_cast(-1.0); + PhasorDynamics::Converter::Reecb bad_reecb_model(&bus, bad_reecb); + success *= (bad_reecb_model.verify() > 0); + + return success.report(__func__); + } + + TestOutcome reecbSignalsInitializationAndResidual() + { + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + using Ext = PhasorDynamics::Converter::ReecbExternalVariables; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + ScalarT pe_value{0.6}; + ScalarT qgen_value{0.2}; + ScalarT iqcmd_value{0.2}; + ScalarT ipcmd_value{0.6}; + IdxT pe_index = 20; + IdxT qgen_index = 21; + IdxT iqcmd_index = 22; + IdxT ipcmd_index = 23; + + PhasorDynamics::SignalNode pe_node; + PhasorDynamics::SignalNode qgen_node; + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + pe_node.set(&pe_value, &pe_index); + qgen_node.set(&qgen_value, &qgen_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + auto data = makeReecbData(); + data.parameters[PhasorDynamics::Converter::ReecbParameters::QFlag] = static_cast(1); + PhasorDynamics::Converter::Reecb reecb(&bus, data); + reecb.getSignals().template attachSignalNode(&pe_node); + reecb.getSignals().template attachSignalNode(&qgen_node); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + iqcmd_node.init(qgen_value); + ipcmd_node.init(pe_value); + success *= (reecb.initialize() == 0); + success *= (reecb.tagDifferentiable() == 0); + success *= (reecb.evaluateResidual() == 0); + + success *= isEqual(reecb.y()[index(Var::VMEAS)], static_cast(1.0), kTol); + success *= isEqual(reecb.y()[index(Var::PMEAS)], pe_value, kTol); + success *= isEqual(reecb.y()[index(Var::QREF)], qgen_value, kTol); + success *= isEqual(iqcmd_node.read(), reecb.y()[index(Var::IQCMD)], kTol); + success *= isEqual(ipcmd_node.read(), reecb.y()[index(Var::IPCMD)], kTol); + success *= (reecb.tag()[index(Var::VMEAS)] == false); + success *= (reecb.tag()[index(Var::PMEAS)] == true); + + for (size_t i = 0; i < reecb.getResidual().size(); ++i) + { + success *= isEqual(reecb.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(reecb.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome rejectsHalfConnectedElectricalFeedback() + { + using Ext = PhasorDynamics::Converter::ReecbExternalVariables; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + ScalarT signal_value{0.6}; + IdxT signal_index = 24; + + PhasorDynamics::SignalNode signal_node; + signal_node.set(&signal_value, &signal_index); + + { + PhasorDynamics::Converter::Reecb reecb(&bus, makeReecbData()); + reecb.getSignals().template attachSignalNode(&signal_node); + success *= (reecb.verify() > 0); + } + + { + PhasorDynamics::Converter::Reecb reecb(&bus, makeReecbData()); + reecb.getSignals().template attachSignalNode(&signal_node); + success *= (reecb.verify() > 0); + } + + return success.report(__func__); + } + + TestOutcome reecbCommandSignalInitialization() + { + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + ScalarT iqcmd_value{0.2}; + ScalarT ipcmd_value{0.6}; + IdxT iqcmd_index = 22; + IdxT ipcmd_index = 23; + + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + PhasorDynamics::Converter::Reecb reecb(&bus, makeReecbData()); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + iqcmd_node.init(static_cast(0.2)); + ipcmd_node.init(static_cast(0.6)); + success *= (reecb.initialize() == 0); + success *= (reecb.tagDifferentiable() == 0); + success *= (reecb.evaluateResidual() == 0); + + success *= isEqual(reecb.y()[index(Var::PMEAS)], static_cast(0.6), kTol); + success *= isEqual(reecb.y()[index(Var::QREF)], static_cast(0.2), kTol); + success *= isEqual(reecb.y()[index(Var::PORD)], static_cast(0.6), kTol); + success *= isEqual(reecb.y()[index(Var::IPCMD)], static_cast(0.6), kTol); + success *= isEqual(reecb.y()[index(Var::IQCMD)], static_cast(0.2), kTol); + + for (size_t i = 0; i < reecb.getResidual().size(); ++i) + { + success *= isEqual(reecb.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(reecb.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome reecbElectricalFeedbackUsesMvaBase() + { + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + using Params = PhasorDynamics::Converter::ReecbParameters; + using Ext = PhasorDynamics::Converter::ReecbExternalVariables; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeReecbData(); + data.parameters[Params::mva] = static_cast(50.0); + + ScalarT pe_value{0.25}; + ScalarT qgen_value{0.05}; + ScalarT iqcmd_value{0.1}; + ScalarT ipcmd_value{0.5}; + IdxT pe_index = 25; + IdxT qgen_index = 26; + IdxT iqcmd_index = 27; + IdxT ipcmd_index = 28; + + PhasorDynamics::SignalNode pe_node; + PhasorDynamics::SignalNode qgen_node; + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + pe_node.set(&pe_value, &pe_index); + qgen_node.set(&qgen_value, &qgen_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + PhasorDynamics::Converter::Reecb reecb(&bus, data); + reecb.getSignals().template attachSignalNode(&pe_node); + reecb.getSignals().template attachSignalNode(&qgen_node); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + iqcmd_node.init(static_cast(0.1)); + ipcmd_node.init(static_cast(0.5)); + success *= (reecb.initialize() == 0); + success *= (reecb.evaluateResidual() == 0); + + success *= isEqual(reecb.y()[index(Var::PMEAS)], static_cast(0.5), kTol); + success *= isEqual(reecb.y()[index(Var::QREF)], static_cast(0.1), kTol); + success *= isEqual(reecb.y()[index(Var::PORD)], static_cast(0.5), kTol); + success *= isEqual(reecb.y()[index(Var::IPCMD)], static_cast(0.5), kTol); + success *= isEqual(reecb.y()[index(Var::IQCMD)], static_cast(0.1), kTol); + + return success.report(__func__); + } + + TestOutcome reecbReferenceFallbackAtAngle() + { + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + using Params = PhasorDynamics::Converter::ReecbParameters; + + TestStatus success = true; + + PhasorDynamics::Bus bus(0.8, 0.6); + bus.allocate(); + bus.initialize(); + + auto data = makeReecbData(); + data.parameters[Params::PfFlag] = static_cast(1); + data.parameters[Params::Qmin] = static_cast(-2.0); + data.parameters[Params::Qmax] = static_cast(2.0); + data.parameters[Params::Pmax] = static_cast(2.0); + data.parameters[Params::Imax] = static_cast(2.0); + + ScalarT iqcmd_value{-0.2}; + ScalarT ipcmd_value{1.6}; + IdxT iqcmd_index = 22; + IdxT ipcmd_index = 23; + + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + PhasorDynamics::Converter::Reecb reecb(&bus, data); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + iqcmd_node.init(static_cast(-0.2)); + ipcmd_node.init(static_cast(1.6)); + success *= (reecb.initialize() == 0); + success *= (reecb.tagDifferentiable() == 0); + success *= (reecb.evaluateResidual() == 0); + + success *= isEqual(reecb.y()[index(Var::PORD)], static_cast(1.6), kTol); + success *= isEqual(reecb.y()[index(Var::QREF)], static_cast(-0.2), kTol); + success *= isEqual(reecb.y()[index(Var::IPCMD)], static_cast(1.6), kTol); + success *= isEqual(reecb.y()[index(Var::IQCMD)], static_cast(-0.2), kTol); + + for (size_t i = 0; i < reecb.getResidual().size(); ++i) + { + success *= isEqual(reecb.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(reecb.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome initializesSaturatedStartConsistently() + { + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + using Params = PhasorDynamics::Converter::ReecbParameters; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.13, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeReecbData(); + data.parameters[Params::QFlag] = static_cast(0); + data.parameters[Params::Pqflag] = static_cast(0); + data.parameters[Params::Vmin] = static_cast(0.9); + data.parameters[Params::Vmax] = static_cast(1.05); + data.parameters[Params::Kvp] = static_cast(10.0); + data.parameters[Params::Kvi] = static_cast(60.0); + data.parameters[Params::Vup] = static_cast(99.0); + data.parameters[Params::Vdip] = static_cast(-99.0); + data.parameters[Params::Imax] = static_cast(1.1); + + ScalarT iqcmd_value{0.15 / 1.13}; + ScalarT ipcmd_value{0.5 / 1.13}; + IdxT iqcmd_index = 22; + IdxT ipcmd_index = 23; + + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + PhasorDynamics::Converter::Reecb reecb(&bus, data); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + iqcmd_node.init(static_cast(0.15 / 1.13)); + ipcmd_node.init(static_cast(0.5 / 1.13)); + success *= (reecb.initialize() == 0); + + const ScalarT piv_arg = static_cast(10.0) * reecb.y()[index(Var::EPIV)] + + reecb.y()[index(Var::XPIV)]; + success *= (piv_arg < -reecb.y()[index(Var::IQMAX)]); + + success *= (reecb.evaluateResidual() == 0); + for (size_t i = 0; i < reecb.getResidual().size(); ++i) + { + success *= isEqual(reecb.getResidual()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + TestOutcome reecbSaturatedVoltagePiJacobianFinite() + { + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + using Params = PhasorDynamics::Converter::ReecbParameters; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.13, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeReecbData(); + data.parameters[Params::QFlag] = static_cast(0); + data.parameters[Params::Pqflag] = static_cast(0); + data.parameters[Params::Vmin] = static_cast(0.9); + data.parameters[Params::Vmax] = static_cast(1.05); + data.parameters[Params::Kvp] = static_cast(10.0); + data.parameters[Params::Kvi] = static_cast(60.0); + data.parameters[Params::Vup] = static_cast(99.0); + data.parameters[Params::Vdip] = static_cast(-99.0); + data.parameters[Params::Imax] = static_cast(1.1); + + ScalarT iqcmd_value{0.15 / 1.13}; + ScalarT ipcmd_value{0.5 / 1.13}; + IdxT iqcmd_index = 22; + IdxT ipcmd_index = 23; + + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + PhasorDynamics::Converter::Reecb reecb(&bus, data); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + iqcmd_node.init(static_cast(0.15 / 1.13)); + ipcmd_node.init(static_cast(0.5 / 1.13)); + success *= (reecb.initialize() == 0); + success *= (reecb.tagDifferentiable() == 0); + success *= (reecb.evaluateResidual() == 0); + success *= (reecb.evaluateJacobian() == 0); + + const ScalarT piv_arg = static_cast(10.0) * reecb.y()[index(Var::EPIV)] + + reecb.y()[index(Var::XPIV)]; + success *= (piv_arg < -reecb.y()[index(Var::IQMAX)]); + + auto jacobian_entries = reecb.getJacobian().getEntries(false); + const auto& values = std::get<2>(jacobian_entries); + success *= (!values.empty()); + for (const auto value : values) + { + success *= std::isfinite(value); + } + + return success.report(__func__); + } +#endif + + TestOutcome jsonParseAndSystemAssembly() + { + TestStatus success = true; + + std::istringstream input(R"json( +{ + "header": { + "format_version": 0, + "format_revision": 1, + "case_name": "renewable electrical control", + "case_description": "REECB parser test", + "case_comments": "", + "freq_base": 60.0, + "va_base": 100000000.0 + }, + "buses": [ + { "number": 1, "class": "bus", "name": "Bus 1", "init": { "Vr": 1.0, "Vi": 0.0 }, "v_base": 1.0 } + ], + "signals": [ + { "signal_id": 12, "name": "Iqcmd" }, + { "signal_id": 13, "name": "Ipcmd" } + ], + "devices": [ + { + "class": "Reecb", + "ports": { "bus": 1, "iqcmd": 12, "ipcmd": 13 }, + "id": "REE1", + "params": { + "mva": 100.0, "PfFlag": 0, "VFlag": 1, "QFlag": 0, "Pqflag": 1, + "Trv": 0.0, "Tp": 0.02, "Vdip": 0.7, "Vup": 1.2, + "dbd1": -0.01, "dbd2": 0.01, "kqv": 0.0, "Iql1": -1.0, "Iqh1": 1.0, + "Qmax": 1.0, "Qmin": -1.0, "Kqp": 1.0, "Kqi": 0.0, + "Vmax": 1.2, "Vmin": 0.8, "Kvp": 1.0, "Kvi": 0.0, + "Tiq": 0.02, "Tpord": 0.02, "dPmax": 1.0, "dPmin": -1.0, + "Pmax": 1.0, "Pmin": 0.0, "Imax": 2.0 + } + } + ] +} +)json"); + + auto data = PhasorDynamics::parseSystemModelData(input); + success *= (data.reecb.size() == 1); + success *= (std::get( + data.reecb[0].parameters.at(PhasorDynamics::Converter::ReecbParameters::Pqflag)) + == static_cast(1)); + success *= (std::get( + data.reecb[0].parameters.at(PhasorDynamics::Converter::ReecbParameters::mva)) + == static_cast(100.0)); + + PhasorDynamics::SystemModel system(data); + success *= (system.allocate() == 0); + success *= (system.initialize() == 0); + success *= (system.evaluateResidual() == 0); + success *= (system.size() == 27); + + return success.report(__func__); + } + + private: + static size_t index(PhasorDynamics::Converter::ReecbInternalVariables variable) + { + return static_cast(variable); + } + + auto makeReecbData() -> PhasorDynamics::Converter::ReecbData + { + using Params = PhasorDynamics::Converter::ReecbParameters; + using Mon = PhasorDynamics::Converter::ReecbMonitorableVariables; + + PhasorDynamics::Converter::ReecbData data; + data.device_class = "Reecb"; + data.disambiguation_string = "reecb_test"; + data.monitored_variables.insert(Mon::iqcmd); + data.monitored_variables.insert(Mon::ipcmd); + + data.parameters[Params::mva] = static_cast(100.0); + data.parameters[Params::PfFlag] = static_cast(0); + data.parameters[Params::VFlag] = static_cast(1); + data.parameters[Params::QFlag] = static_cast(0); + data.parameters[Params::Pqflag] = static_cast(1); + data.parameters[Params::Trv] = static_cast(0.0); + data.parameters[Params::Tp] = static_cast(0.02); + data.parameters[Params::Vdip] = static_cast(0.7); + data.parameters[Params::Vup] = static_cast(1.2); + data.parameters[Params::dbd1] = static_cast(-0.01); + data.parameters[Params::dbd2] = static_cast(0.01); + data.parameters[Params::kqv] = static_cast(0.0); + data.parameters[Params::Iql1] = static_cast(-1.0); + data.parameters[Params::Iqh1] = static_cast(1.0); + data.parameters[Params::Qmax] = static_cast(1.0); + data.parameters[Params::Qmin] = static_cast(-1.0); + data.parameters[Params::Kqp] = static_cast(1.0); + data.parameters[Params::Kqi] = static_cast(0.0); + data.parameters[Params::Vmax] = static_cast(1.2); + data.parameters[Params::Vmin] = static_cast(0.8); + data.parameters[Params::Kvp] = static_cast(1.0); + data.parameters[Params::Kvi] = static_cast(0.0); + data.parameters[Params::Tiq] = static_cast(0.02); + data.parameters[Params::Tpord] = static_cast(0.02); + data.parameters[Params::dPmax] = static_cast(1.0); + data.parameters[Params::dPmin] = static_cast(-1.0); + data.parameters[Params::Pmax] = static_cast(1.0); + data.parameters[Params::Pmin] = static_cast(0.0); + data.parameters[Params::Imax] = static_cast(2.0); + + return data; + } + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runConverterReecbTests.cpp b/tests/UnitTests/PhasorDynamics/runConverterReecbTests.cpp new file mode 100644 index 000000000..0645f6133 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runConverterReecbTests.cpp @@ -0,0 +1,22 @@ +#include "ConverterReecbTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::ConverterReecbTests test; + + result += test.constructionAndValidation(); + result += test.reecbSignalsInitializationAndResidual(); + result += test.rejectsHalfConnectedElectricalFeedback(); + result += test.reecbCommandSignalInitialization(); + result += test.reecbElectricalFeedbackUsesMvaBase(); + result += test.reecbReferenceFallbackAtAngle(); + result += test.initializesSaturatedStartConsistently(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.reecbSaturatedVoltagePiJacobianFinite(); +#endif + result += test.jsonParseAndSystemAssembly(); + + return result.summary(); +}