From 44ca9b30d9fd00a6821fdf932f46c1a2eeab5dab Mon Sep 17 00:00:00 2001 From: lukelowry Date: Thu, 28 May 2026 19:01:57 -0500 Subject: [PATCH] Add HYGOV phasor dynamics governor --- CHANGELOG.md | 1 + GridKit/CommonMath.hpp | 30 +- GridKit/CommonMath.md | 16 +- .../Model/PhasorDynamics/ComponentLibrary.hpp | 1 + .../PhasorDynamics/Governor/CMakeLists.txt | 1 + .../Governor/HYGOV/CMakeLists.txt | 54 ++ .../PhasorDynamics/Governor/HYGOV/Hygov.cpp | 27 + .../PhasorDynamics/Governor/HYGOV/Hygov.hpp | 154 +++++ .../Governor/HYGOV/HygovData.hpp | 91 +++ .../HYGOV/HygovDependencyTracking.cpp | 27 + .../Governor/HYGOV/HygovEnzyme.cpp | 111 ++++ .../Governor/HYGOV/HygovImpl.hpp | 553 ++++++++++++++++++ .../PhasorDynamics/Governor/HYGOV/README.md | 291 +++++++++ .../Model/PhasorDynamics/Governor/README.md | 1 + GridKit/Model/PhasorDynamics/INPUT_FORMAT.md | 1 + GridKit/Model/PhasorDynamics/SystemModel.hpp | 32 + .../Model/PhasorDynamics/SystemModelData.hpp | 3 + .../SystemModelDataJSONParser.hpp | 6 + docs/Figures/PhasorDynamics/HYGOV_diagram.png | Bin 0 -> 72381 bytes .../Math/SmoothnessIndicatorTests.hpp | 50 +- .../Math/runSmoothnessIndicatorTests.cpp | 3 +- tests/UnitTests/PhasorDynamics/CMakeLists.txt | 10 + .../PhasorDynamics/GovernorHygovTests.hpp | 363 ++++++++++++ .../PhasorDynamics/runGovernorHygovTests.cpp | 16 + 24 files changed, 1822 insertions(+), 20 deletions(-) create mode 100644 GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt create mode 100644 GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp create mode 100644 GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp create mode 100644 GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp create mode 100644 GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp create mode 100644 GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp create mode 100644 GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp create mode 100644 GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md create mode 100755 docs/Figures/PhasorDynamics/HYGOV_diagram.png create mode 100644 tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp create mode 100644 tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index 066786e09..2eaacb7b3 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 portable Vector class to GridKit +- Added `HYGOV` governor model implementation for PhasorDynamics. ## v0.1 diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index 9bb3a5c3c..daaa11e94 100644 --- a/GridKit/CommonMath.hpp +++ b/GridKit/CommonMath.hpp @@ -153,7 +153,31 @@ namespace GridKit } /** - * @brief Smooth two-sided deadband function + * @brief Smooth Type 1 no-offset two-sided deadband function + * + * Smooth approximation to a deadband that returns zero inside the band and + * passes the input through unchanged outside the band. + * + * @tparam ScalarT - scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] x - Input signal + * @param[in] lower - Lower breakpoint + * @param[in] upper - Upper breakpoint + * @return Smooth no-offset deadbanded value + */ + template + __attribute__((always_inline)) inline ScalarT deadband1( + const ScalarT x, + const RealT lower, + const RealT upper) + { + assert(lower <= upper); + return x * (sigmoid(lower - x) + sigmoid(x - upper)); + } + + /** + * @brief Smooth Type 2 offset two-sided deadband function * * Smooth approximation to x - min(max(x, lower), upper), composed from the * smooth ramp function. @@ -164,10 +188,10 @@ namespace GridKit * @param[in] x - Input signal * @param[in] lower - Lower breakpoint * @param[in] upper - Upper breakpoint - * @return Smooth deadbanded value + * @return Smooth offset deadbanded value */ template - __attribute__((always_inline)) inline ScalarT deadband( + __attribute__((always_inline)) inline ScalarT deadband2( const ScalarT x, const RealT lower, const RealT upper) diff --git a/GridKit/CommonMath.md b/GridKit/CommonMath.md index a654487a2..29a9a7901 100644 --- a/GridKit/CommonMath.md +++ b/GridKit/CommonMath.md @@ -49,7 +49,8 @@ The scale $\mu=4\cdot f_{\text{sync}}=240$ is chosen so $\sigma$ behaves like a | `max` | Smooth binary maximum | `REECA` | | `min` | Smooth binary minimum | `REECA` | | `clamp` | Bounded saturation | `IEEEST` | -| `deadband` | Signed two-sided deadband | `REECA` | +| `deadband1` | Type 1 no-offset signed two-sided deadband | - | +| `deadband2` | Type 2 offset signed two-sided deadband | `REECA` | | `slew` | Symmetric slew-rate limiter | - | | `linseg` | Saturated linear segment contribution | `REGCA`, `REECA` | | `above` | Above-lower-limit indicator | - | @@ -84,7 +85,15 @@ x & \ell \le x \le u \\ u & x > u \end{cases} \\ -\text{deadband}(x;\ell,u) +\text{deadband1}(x;\ell,u) +&= +\begin{cases} +x & x < \ell \\ +0 & \ell \le x \le u \\ +x & x > u +\end{cases} +\\ +\text{deadband2}(x;\ell,u) &= \begin{cases} x-\ell & x < \ell \\ @@ -156,7 +165,8 @@ f & x \ge u \land f < 0 \\ \text{max}(x,y) &= y+\rho(x-y) \\ \text{min}(x,y) &= x-\rho(x-y) \\ \text{clamp}(x;\ell,u) &= \ell+\rho(x-\ell)-\rho(x-u) \\ -\text{deadband}(x;\ell,u) &= \rho(x-u)-\rho(\ell-x) \\ +\text{deadband1}(x;\ell,u) &= x\left[\sigma(\ell-x)+\sigma(x-u)\right] \\ +\text{deadband2}(x;\ell,u) &= \rho(x-u)-\rho(\ell-x) \\ \text{slew}(f;r) &= -r+\rho(f+r)-\rho(f-r) \\ \text{linseg}(x;a,b,h) &= \dfrac{h}{b-a}\left[\rho(x-a)-\rho(x-b)\right] \\ \text{above}(x;\ell) &= \sigma(x-\ell) \\ diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index e48a8a9c5..0a033850a 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt index 7c7269784..fbe71c740 100644 --- a/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt @@ -4,3 +4,4 @@ # ]] add_subdirectory(Tgov1) +add_subdirectory(HYGOV) diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt new file mode 100644 index 000000000..1719101b0 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Hygov.hpp HygovData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_governor_hygov + SOURCES HygovEnzyme.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_governor_hygov + SOURCES Hygov.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_governor_hygov_dependency_tracking + SOURCES HygovDependencyTracking.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_governor_hygov) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_governor_hygov_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp new file mode 100644 index 000000000..bdb7e15c0 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp @@ -0,0 +1,27 @@ +/** + * @file Hygov.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the HYGOV governor model. + */ + +#include "HygovImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Hygov; + template class Hygov; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp new file mode 100644 index 000000000..19adabede --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp @@ -0,0 +1,154 @@ +/** + * @file Hygov.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the HYGOV governor model. + */ + +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class SignalNode; + + namespace Governor + { + /// Internal variables of a `Hygov`. + enum class HygovInternalVariables : size_t + { + XN, ///< Speed lead-lag denominator state + XF, ///< Governor error filter output + C, ///< Desired-gate position + G, ///< Gate position + Q, ///< Turbine flow + OMEGADB, ///< Deadbanded speed deviation + YOMEGA, ///< Lead-lag-conditioned speed signal + EF, ///< Governor error into the filter + FC, ///< Desired-gate derivative target + RC, ///< Rate-limited desired-gate derivative target + PGV, ///< Nonlinear gate-to-power curve output + PGVSAFE, ///< Safe gate-to-power value + H, ///< Turbine head + PMECH, ///< Mechanical-power output + MAXIMUM, + }; + + /// External variables of a `Hygov`. + enum class HygovExternalVariables : size_t + { + OMEGA, ///< Machine speed deviation + PREF, ///< Active-power/load reference + PAUX, ///< Auxiliary power input + MAXIMUM, + }; + + template + class Hygov : public Component + { + using Component::alpha_; + using Component::f_; + using Component::gridkit_component_id_; + using Component::J_; + using Component::J_cols_buffer_; + using Component::J_rows_buffer_; + using Component::J_vals_buffer_; + using Component::residual_indices_; + using Component::size_; + using Component::tag_; + using Component::variable_indices_; + using Component::wb_; + using Component::y_; + using Component::yp_; + + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename Component::RealT; + using SignalT = SignalNode; + using ModelDataT = HygovData; + using MonitorT = Model::VariableMonitor; + + Hygov(); + Hygov(const ModelDataT& data); + ~Hygov(); + + 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 ModelDataT& data); + void initializeMonitor(); + + ScalarT gatePower(ScalarT gate) const; + RealT invertGatePower(RealT pgv) const; + ScalarT toComponentBase(ScalarT value) const; + ScalarT toPmechBase(ScalarT value) const; + + RealT Trate_{0}; + RealT pmech_mva_base_{0}; + RealT Rperm_{0}; + RealT Rtemp_{0}; + RealT Tr_{0}; + RealT Tf_{0}; + RealT Tg_{0}; + RealT Velm_{0}; + RealT Gmax_{0}; + RealT Gmin_{0}; + RealT Tw_{0}; + RealT At_{0}; + RealT Dturb_{0}; + RealT Qnl_{0}; + RealT Tn_{0}; + RealT Tnp_{0}; + RealT leadlag_gain_{0}; + RealT db1_{0}; + RealT db2_{0}; + RealT Hdam_{1}; + std::array Gv_{}; + std::array Pgv_{}; + + IdxT parameter_error_count_{0}; + + ScalarT pref_set_{0}; + ScalarT paux_set_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp new file mode 100644 index 000000000..9e2ebfda9 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp @@ -0,0 +1,91 @@ +/** + * @file HygovData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the HYGOV governor model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + /// Parameter keys for the HYGOV governor model. + enum class HygovParameters + { + Trate, ///< Turbine-rating power base + mva, ///< MVA base of the pmech signal + Rperm, ///< Permanent droop + Rtemp, ///< Temporary droop + Tr, ///< Temporary-droop reset time constant + Tf, ///< Governor error filter time constant + Tg, ///< Gate servo time constant + Velm, ///< Maximum desired-gate velocity magnitude + Gmax, ///< Maximum desired-gate position + Gmin, ///< Minimum desired-gate position + Tw, ///< Water inertia time constant + At, ///< Turbine gain + Dturb, ///< Turbine damping coefficient + Qnl, ///< No-load flow at nominal head + Tn, ///< Speed lead-lag numerator time constant + Tnp, ///< Speed lead-lag denominator time constant + db1, ///< Type 1 speed deadband threshold + db2, ///< Unsupported mechanical backlash deadband + Hdam, ///< Head available at dam + Gv0, ///< Gate point 0 + Gv1, ///< Gate point 1 + Gv2, ///< Gate point 2 + Gv3, ///< Gate point 3 + Gv4, ///< Gate point 4 + Gv5, ///< Gate point 5 + Pgv0, ///< Power point 0 + Pgv1, ///< Power point 1 + Pgv2, ///< Power point 2 + Pgv3, ///< Power point 3 + Pgv4, ///< Power point 4 + Pgv5 ///< Power point 5 + }; + + /// Ports for the HYGOV governor model. + enum class HygovPorts + { + bus, ///< Optional terminal bus ID for case-format compatibility + speed, ///< Machine speed-deviation signal ID + pmech, ///< Mechanical-power output signal ID + pref, ///< Optional active-power/load reference signal ID + paux ///< Optional auxiliary power input signal ID + }; + + /// Variables available through the monitor interface. + enum class HygovMonitorableVariables + { + pmech, ///< Mechanical power output + leadlag, ///< Speed lead-lag denominator state + filter, ///< Governor error filter output + desiredgate, ///< Desired-gate position + gate, ///< Gate position + flow, ///< Turbine flow + head, ///< Turbine head + pgv ///< Nonlinear gate-to-power curve output + }; + + template + struct HygovData : public ComponentData + { + HygovData() = default; + + using Parameters = HygovParameters; + using Ports = HygovPorts; + using MonitorableVariables = HygovMonitorableVariables; + }; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp new file mode 100644 index 000000000..760bac957 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file HygovDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the HYGOV governor model. + */ + +#include "HygovImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Hygov; + template class Hygov; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp new file mode 100644 index 000000000..7c5ff0aac --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp @@ -0,0 +1,111 @@ +/** + * @file HygovEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the HYGOV governor model. + */ + +#include + +#include "HygovImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + namespace + { + template + real_type deadband1Derivative(real_type x, real_type lower, real_type upper) + { + using RealT = real_type; + static constexpr RealT MU = static_cast(240.0); + + const RealT lower_sig = Math::sigmoid(lower - x); + const RealT upper_sig = Math::sigmoid(x - upper); + return lower_sig + upper_sig + + x * MU + * (upper_sig * (ONE - upper_sig) + - lower_sig * (ONE - lower_sig)); + } + } // namespace + + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << 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::Governor::Hygov; + 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_); + + const auto OMEGADB = static_cast(HygovInternalVariables::OMEGADB); + const auto EF = static_cast(HygovInternalVariables::EF); + const auto PMECH = static_cast(HygovInternalVariables::PMECH); + const auto G = static_cast(HygovInternalVariables::G); + + const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); + const auto PREF = static_cast(HygovExternalVariables::PREF); + const auto PAUX = static_cast(HygovExternalVariables::PAUX); + + IdxT nnz = 0; + auto storeSignalDerivative = [&](size_t row, size_t col, RealT value) + { + if (ws_indices_[col] == INVALID_INDEX) + { + return; + } + + J_rows_buffer_[nnz] = residual_indices_.at(row); + J_cols_buffer_[nnz] = ws_indices_[col]; + J_vals_buffer_[nnz] = value; + ++nnz; + }; + + storeSignalDerivative(OMEGADB, + OMEGA, + deadband1Derivative(static_cast(ws_[OMEGA]), + -db1_, + db1_)); + storeSignalDerivative(EF, PREF, ONE); + storeSignalDerivative(EF, PAUX, ONE); + storeSignalDerivative(PMECH, + OMEGA, + -Dturb_ * static_cast(y_[G]) * Trate_ / pmech_mva_base_); + + J_.setValues(1.0, J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, nnz); + return 0; + } + + template class Hygov; + template class Hygov; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp new file mode 100644 index 000000000..a24158e3d --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp @@ -0,0 +1,553 @@ +/** + * @file HygovImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the HYGOV governor model. + */ + +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + using Log = ::GridKit::Utilities::Logger; + + template + Hygov::Hygov() + { + size_ = static_cast(HygovInternalVariables::MAXIMUM); + } + + template + Hygov::Hygov(const ModelDataT& data) + : monitor_(std::make_unique(data)) + { + initModelParams(data); + initializeMonitor(); + size_ = static_cast(HygovInternalVariables::MAXIMUM); + } + + template + Hygov::~Hygov() + { + } + + template + scalar_type Hygov::toComponentBase(scalar_type value) const + { + return value * pmech_mva_base_ / Trate_; + } + + template + scalar_type Hygov::toPmechBase(scalar_type value) const + { + return value * Trate_ / pmech_mva_base_; + } + + template + void Hygov::initModelParams(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + + parameter_error_count_ = 0; + + auto load_required_real = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Hygov: 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() << "Hygov: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + auto load_optional_real = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + 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() << "Hygov: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + load_required_real(Params::Trate, Trate_, "Trate"); + load_required_real(Params::mva, pmech_mva_base_, "mva"); + load_required_real(Params::Rperm, Rperm_, "Rperm"); + load_required_real(Params::Rtemp, Rtemp_, "Rtemp"); + load_required_real(Params::Tr, Tr_, "Tr"); + load_required_real(Params::Tf, Tf_, "Tf"); + load_required_real(Params::Tg, Tg_, "Tg"); + load_required_real(Params::Velm, Velm_, "Velm"); + load_required_real(Params::Gmax, Gmax_, "Gmax"); + load_required_real(Params::Gmin, Gmin_, "Gmin"); + load_required_real(Params::Tw, Tw_, "Tw"); + load_required_real(Params::At, At_, "At"); + load_required_real(Params::Dturb, Dturb_, "Dturb"); + load_required_real(Params::Qnl, Qnl_, "Qnl"); + load_required_real(Params::Tn, Tn_, "Tn"); + load_required_real(Params::Tnp, Tnp_, "Tnp"); + load_required_real(Params::db1, db1_, "db1"); + load_optional_real(Params::db2, db2_, "db2"); + load_required_real(Params::Hdam, Hdam_, "Hdam"); + + load_required_real(Params::Gv0, Gv_[0], "Gv0"); + load_required_real(Params::Gv1, Gv_[1], "Gv1"); + load_required_real(Params::Gv2, Gv_[2], "Gv2"); + load_required_real(Params::Gv3, Gv_[3], "Gv3"); + load_required_real(Params::Gv4, Gv_[4], "Gv4"); + load_required_real(Params::Gv5, Gv_[5], "Gv5"); + load_required_real(Params::Pgv0, Pgv_[0], "Pgv0"); + load_required_real(Params::Pgv1, Pgv_[1], "Pgv1"); + load_required_real(Params::Pgv2, Pgv_[2], "Pgv2"); + load_required_real(Params::Pgv3, Pgv_[3], "Pgv3"); + load_required_real(Params::Pgv4, Pgv_[4], "Pgv4"); + load_required_real(Params::Pgv5, Pgv_[5], "Pgv5"); + + const bool source_default_curve = + std::all_of(Gv_.begin(), Gv_.end(), [](RealT value) + { return value == ZERO; }) + && std::all_of(Pgv_.begin(), Pgv_.end(), [](RealT value) + { return value == ZERO; }); + if (source_default_curve) + { + Gv_ = {ZERO, + static_cast(0.2), + static_cast(0.4), + static_cast(0.6), + static_cast(0.8), + ONE}; + Pgv_ = Gv_; + } + + leadlag_gain_ = ZERO; + if (Tnp_ > ZERO) + { + leadlag_gain_ = Tn_ / Tnp_; + } + } + + template + const Model::VariableMonitorBase* Hygov::getMonitor() const + { + return monitor_.get(); + } + + template + void Hygov::initializeMonitor() + { + using Variable = typename ModelDataT::MonitorableVariables; + auto index = [](HygovInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::pmech, [this, index] + { return y_[index(HygovInternalVariables::PMECH)]; }); + monitor_->set(Variable::leadlag, [this, index] + { return y_[index(HygovInternalVariables::XN)]; }); + monitor_->set(Variable::filter, [this, index] + { return y_[index(HygovInternalVariables::XF)]; }); + monitor_->set(Variable::desiredgate, [this, index] + { return y_[index(HygovInternalVariables::C)]; }); + monitor_->set(Variable::gate, [this, index] + { return y_[index(HygovInternalVariables::G)]; }); + monitor_->set(Variable::flow, [this, index] + { return y_[index(HygovInternalVariables::Q)]; }); + monitor_->set(Variable::head, [this, index] + { return y_[index(HygovInternalVariables::H)]; }); + monitor_->set(Variable::pgv, [this, index] + { return y_[index(HygovInternalVariables::PGV)]; }); + } + + template + int Hygov::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Hygov::allocate() + { + size_ = static_cast(HygovInternalVariables::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_.clear(); + + auto signal_size = static_cast(HygovExternalVariables::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(HygovInternalVariables::PMECH)], + &(this->getVariableIndex(static_cast(HygovInternalVariables::PMECH)))); + } + + return 0; + } + + template + int Hygov::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Hygov: " << message << '\n'; + ret += 1; + } + }; + + check(Trate_ > ZERO, "Trate must be positive"); + check(pmech_mva_base_ > ZERO, "mva must be positive"); + check(Rperm_ > ZERO, "Rperm must be positive"); + check(Rtemp_ > ZERO, "Rtemp must be positive"); + check(Tr_ > ZERO, "Tr must be positive"); + check(Tf_ > ZERO, "Tf must be positive"); + check(Tg_ >= ZERO, "Tg must be non-negative"); + check(Tw_ > ZERO, "Tw must be positive"); + check(Tn_ >= ZERO, "Tn must be non-negative"); + check(Tnp_ >= ZERO, "Tnp must be non-negative"); + check(Tnp_ > ZERO || Tn_ == ZERO, "Tn must be zero when Tnp is zero"); + check(Velm_ > ZERO, "Velm must be positive"); + check(Gmin_ <= Gmax_, "Gmin must be less than or equal to Gmax"); + check(At_ > ZERO, "At must be positive"); + check(Dturb_ >= ZERO, "Dturb must be non-negative"); + check(Qnl_ >= ZERO, "Qnl must be non-negative"); + check(db1_ >= ZERO, "db1 must be non-negative"); + check(db2_ == ZERO, "nonzero db2 is not supported"); + check(Hdam_ > ZERO, "Hdam must be positive"); + + for (size_t i = 1; i < Gv_.size(); ++i) + { + check(Gv_[i - 1] < Gv_[i], "Gv points must be strictly increasing"); + check(Pgv_[i - 1] <= Pgv_[i], "Pgv points must be non-decreasing"); + } + check(Pgv_[0] >= ZERO, "Pgv0 must be non-negative"); + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Hygov: omega signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Hygov: pref signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Hygov: paux signal attached with no linked source\n"; + ret += 1; + } + + return ret; + } + + template + scalar_type Hygov::gatePower(scalar_type gate) const + { + return ScalarT{Pgv_[0]} + + Math::linseg(gate, Gv_[0], Gv_[1], Pgv_[1] - Pgv_[0]) + + Math::linseg(gate, Gv_[1], Gv_[2], Pgv_[2] - Pgv_[1]) + + Math::linseg(gate, Gv_[2], Gv_[3], Pgv_[3] - Pgv_[2]) + + Math::linseg(gate, Gv_[3], Gv_[4], Pgv_[4] - Pgv_[3]) + + Math::linseg(gate, Gv_[4], Gv_[5], Pgv_[5] - Pgv_[4]); + } + + template + typename Hygov::RealT + Hygov::invertGatePower( + typename Hygov::RealT pgv) const + { + static constexpr RealT tol = static_cast(1.0e-10); + + if (std::abs(pgv - Pgv_[0]) <= tol) + { + return Gv_[0]; + } + + for (size_t i = 0; i < 5; ++i) + { + if (Pgv_[i + 1] <= Pgv_[i]) + { + continue; + } + + if (Pgv_[i] - tol <= pgv && pgv <= Pgv_[i + 1] + tol) + { + const RealT fraction = (pgv - Pgv_[i]) / (Pgv_[i + 1] - Pgv_[i]); + return Gv_[i] + fraction * (Gv_[i + 1] - Gv_[i]); + } + } + + return std::numeric_limits::quiet_NaN(); + } + + template + int Hygov::initialize() + { + if (parameter_error_count_ > 0 || verify() > 0) + { + Log::error() << "Hygov: cannot initialize with invalid configuration\n"; + return 1; + } + + const auto XN = static_cast(HygovInternalVariables::XN); + const auto XF = static_cast(HygovInternalVariables::XF); + const auto C = static_cast(HygovInternalVariables::C); + const auto G = static_cast(HygovInternalVariables::G); + const auto Q = static_cast(HygovInternalVariables::Q); + const auto OMEGADB = static_cast(HygovInternalVariables::OMEGADB); + const auto YOMEGA = static_cast(HygovInternalVariables::YOMEGA); + const auto EF = static_cast(HygovInternalVariables::EF); + const auto FC = static_cast(HygovInternalVariables::FC); + const auto RC = static_cast(HygovInternalVariables::RC); + const auto PGV = static_cast(HygovInternalVariables::PGV); + const auto PGVSAFE = static_cast(HygovInternalVariables::PGVSAFE); + const auto H = static_cast(HygovInternalVariables::H); + const auto PMECH = static_cast(HygovInternalVariables::PMECH); + + ScalarT omega0{ZERO}; + if (signals_.template isAttached()) + { + omega0 = signals_.template readExternalVariable(); + } + + paux_set_ = ScalarT{ZERO}; + if (signals_.template isAttached()) + { + paux_set_ = signals_.template readExternalVariable(); + } + + const ScalarT pmech0 = toComponentBase(y_[PMECH]); + y_[H] = Hdam_; + y_[Q] = Qnl_ + pmech0 / (At_ * y_[H]); + y_[PGV] = y_[Q] / std::sqrt(y_[H]); + y_[PGVSAFE] = Math::max(y_[PGV], static_cast(0.01)); + + const RealT gate0 = invertGatePower(static_cast(y_[PGV])); + if (std::isnan(gate0)) + { + Log::error() << "Hygov: initial Pgv is outside the invertible gate curve\n"; + return 1; + } + + y_[G] = gate0; + y_[C] = y_[G]; + + if (y_[C] < Gmin_ || y_[C] > Gmax_) + { + Log::error() << "Hygov: initialized gate is outside Gmin/Gmax\n"; + return 1; + } + + y_[OMEGADB] = Math::deadband1(omega0, -db1_, db1_); + y_[XN] = y_[OMEGADB]; + y_[YOMEGA] = y_[OMEGADB]; + y_[XF] = ZERO; + y_[EF] = ZERO; + y_[FC] = ZERO; + y_[RC] = ZERO; + y_[PMECH] = toPmechBase(pmech0); + + pref_set_ = y_[EF] - paux_set_ + y_[YOMEGA] + Rperm_ * y_[C]; + if (signals_.template isAttached()) + { + const ScalarT pref0 = + signals_.template readExternalVariable(); + const RealT pref_err = static_cast(pref0 - pref_set_); + if (std::abs(pref_err) > static_cast(1.0e-10)) + { + Log::error() << "Hygov: pref initial condition is not steady state\n"; + return 1; + } + pref_set_ = pref0; + } + + std::fill(yp_.begin(), yp_.end(), ZERO); + return 0; + } + + template + int Hygov::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(HygovInternalVariables::XN)] = (Tnp_ > ZERO); + tag_[static_cast(HygovInternalVariables::XF)] = true; + tag_[static_cast(HygovInternalVariables::C)] = true; + tag_[static_cast(HygovInternalVariables::G)] = (Tg_ > ZERO); + tag_[static_cast(HygovInternalVariables::Q)] = true; + return 0; + } + + template + __attribute__((always_inline)) inline int + Hygov::evaluateInternalResidual( + ScalarT* y, + ScalarT* yp, + [[maybe_unused]] ScalarT* wb, + ScalarT* ws, + ScalarT* f) + { + const auto XN = static_cast(HygovInternalVariables::XN); + const auto XF = static_cast(HygovInternalVariables::XF); + const auto C = static_cast(HygovInternalVariables::C); + const auto G = static_cast(HygovInternalVariables::G); + const auto Q = static_cast(HygovInternalVariables::Q); + const auto OMEGADB = static_cast(HygovInternalVariables::OMEGADB); + const auto YOMEGA = static_cast(HygovInternalVariables::YOMEGA); + const auto EF = static_cast(HygovInternalVariables::EF); + const auto FC = static_cast(HygovInternalVariables::FC); + const auto RC = static_cast(HygovInternalVariables::RC); + const auto PGV = static_cast(HygovInternalVariables::PGV); + const auto PGVSAFE = static_cast(HygovInternalVariables::PGVSAFE); + const auto H = static_cast(HygovInternalVariables::H); + const auto PMECH = static_cast(HygovInternalVariables::PMECH); + + const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); + const auto PREF = static_cast(HygovExternalVariables::PREF); + const auto PAUX = static_cast(HygovExternalVariables::PAUX); + + const ScalarT xn = y[XN]; + const ScalarT xf = y[XF]; + const ScalarT c = y[C]; + const ScalarT g = y[G]; + const ScalarT q = y[Q]; + const ScalarT omegadb = y[OMEGADB]; + const ScalarT yomega = y[YOMEGA]; + const ScalarT ef = y[EF]; + const ScalarT fc = y[FC]; + const ScalarT rc = y[RC]; + const ScalarT pgv = y[PGV]; + const ScalarT pgv_safe = y[PGVSAFE]; + const ScalarT head = y[H]; + const ScalarT pmech = y[PMECH]; + + const ScalarT omega = ws[OMEGA]; + const ScalarT pref = ws[PREF]; + const ScalarT paux = ws[PAUX]; + + f[XN] = -Tnp_ * yp[XN] - xn + omegadb; + f[XF] = -Tf_ * yp[XF] - xf + ef; + f[FC] = -Rtemp_ * Tr_ * fc + xf + Tr_ * yp[XF]; + f[C] = -yp[C] + Math::antiwindup(c, rc, Gmin_, Gmax_); + f[G] = -Tg_ * yp[G] - g + c; + f[Q] = -Tw_ * yp[Q] + Hdam_ - head; + + f[OMEGADB] = -omegadb + Math::deadband1(omega, -db1_, db1_); + f[YOMEGA] = -yomega + xn + leadlag_gain_ * (omegadb - xn); + f[EF] = -ef + pref + paux - yomega - Rperm_ * c; + f[RC] = -rc + Math::clamp(fc, -Velm_, Velm_); + f[PGV] = -pgv + gatePower(g); + f[PGVSAFE] = -pgv_safe + Math::max(pgv, static_cast(0.01)); + f[H] = -head + (q / pgv_safe) * (q / pgv_safe); + const ScalarT pmech_output = At_ * head * (q - Qnl_) - Dturb_ * omega * g; + f[PMECH] = -pmech + toPmechBase(pmech_output); + + return 0; + } + + template + int Hygov::evaluateResidual() + { + const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); + const auto PREF = static_cast(HygovExternalVariables::PREF); + const auto PAUX = static_cast(HygovExternalVariables::PAUX); + + ws_[OMEGA] = ZERO; + ws_[PREF] = pref_set_; + ws_[PAUX] = paux_set_; + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + if (signals_.template isAttached()) + { + ws_[OMEGA] = signals_.template readExternalVariable(); + ws_indices_[OMEGA] = + signals_.template readExternalVariableIndex(); + } + + if (signals_.template isAttached()) + { + ws_[PREF] = signals_.template readExternalVariable(); + ws_indices_[PREF] = + signals_.template readExternalVariableIndex(); + } + + if (signals_.template isAttached()) + { + ws_[PAUX] = signals_.template readExternalVariable(); + ws_indices_[PAUX] = + signals_.template readExternalVariableIndex(); + } + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + return 0; + } + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md b/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md new file mode 100644 index 000000000..d0957ecdd --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md @@ -0,0 +1,291 @@ +# **Hydro Turbine-Governor Model (HYGOV)** + +HYGOV is a hydro turbine-governor model with a temporary-droop governor, gate +servo, and nonlinear single-penstock water-column turbine. In GridKit it is +represented as a governor model that reads machine speed deviation and supplies +mechanical power to the machine. + +Notes: +- HYGOV equations use the `Trate` base; the `pmech` signal uses `mva` and is converted at the signal boundary. +- GridKit requires `Trate > 0`; for PowerWorld `Trate = 0` records, set `Trate` to the connected machine `mva` before running. +- Aside from `db2`, model parameters are required input values; the table values are typical values, not parser defaults. +- HYGOV uses `db1` with CommonMath `deadband1`; HYGOVD-only `dbL`/`dbH` and nonzero `db2` are not modeled. +- PowerWorld fields `Ttur`, `Eps`, Kaplan blade-servo points, `Bmax`, and `Tblade` are not part of this implementation. +- Permanent-droop feedback uses desired gate $c$ (State 2), not gate $g$. + +## Block Diagram + +Standard HYGOV block diagram. + +
+ + + Figure 1: Governor HYGOV model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +
+ +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value | Note +--------------------------------|----------|----------|----------------------------------------------|---------------|------ +$P^{\mathrm{rate}}$ | [MW] | `Trate` | Component-base turbine rating | 100.0 | Required positive value +$S^{\mathrm{pmech}}$ | [MVA] | `mva` | Mechanical-power signal base | 100.0 | Use connected machine `mva` +$R_{\mathrm{perm}}$ | [p.u.] | `Rperm` | Permanent droop | 0.04 | Source diagram label: `R` +$R_{\mathrm{temp}}$ | [p.u.] | `Rtemp` | Temporary droop | 0.3 | Source diagram label: `r` +$T_r$ | [sec] | `Tr` | Temporary-droop reset time constant | 5.0 | +$T_f$ | [sec] | `Tf` | Governor error filter time constant | 0.05 | State 1 +$T_g$ | [sec] | `Tg` | Gate servo time constant | 0.5 | State 3; if zero, $g$ is algebraic +$V_{\mathrm{elm}}$ | [p.u./s] | `Velm` | Maximum desired-gate velocity magnitude | 0.2 | Symmetric rate limit on State 2 +$G^{\max}$ | [p.u.] | `Gmax` | Maximum desired-gate position | 1.0 | +$G^{\min}$ | [p.u.] | `Gmin` | Minimum desired-gate position | 0.0 | +$T_w$ | [sec] | `Tw` | Water inertia time constant | 1.0 | State 4 +$A_t$ | [p.u.] | `At` | Turbine gain | 1.2 | +$D_{\mathrm{turb}}$ | [p.u.] | `Dturb` | Turbine damping coefficient | 0.5 | Multiplied by speed deviation and gate +$q_{\mathrm{NL}}$ | [p.u.] | `Qnl` | No-load flow at nominal head | 0.05 | +$T_n$ | [sec] | `Tn` | Speed lead-lag numerator time constant | 0.0 | Source text labels this "lag" +$T_{np}$ | [sec] | `Tnp` | Speed lead-lag denominator time constant | 0.0 | Source text labels this "lead" +$D_{\omega}$ | [p.u.] | `db1` | Type 1 speed deadband threshold | 0.0 | Source data may provide this in Hz +$H_{\mathrm{dam}}$ | [p.u.] | `Hdam` | Head available at dam | 1.0 | + +The nonlinear gate-to-power curve $N_{\mathrm{GV}}$ is represented by six source points: + +Symbol | Units | JSON | Description | Typical Value | Note +--------------------------------|--------|----------|----------------------------------------------|---------------|------ +$G_V^{(k)}$ | [p.u.] | `Gv0`-`Gv5` | Gate point $k$ of the gain curve | 0.0 | $k=0,\ldots,5$ +$P_{\mathrm{GV}}^{(k)}$ | [p.u.] | `Pgv0`-`Pgv5` | Power point $k$ of the gain curve | 0.0 | $k=0,\ldots,5$ + +If all six `Gv` and all six `Pgv` source values are zero, the source data is treated as +omitting the nonlinear curve. The effective curve used by the equations is the identity +curve with points `(0,0), (0.2,0.2), ..., (1,1)`. + +### Parameter Validation + +Invalid HYGOV parameter sets are rejected by the following checks. If source data preprocessing swaps gate limits, expands gate limits to include the initialized gate, adjusts small positive time constants, converts `db1` from Hz to p.u. speed, or replaces an all-zero source gate curve with the identity curve, apply these checks to the effective values used by the equations. Source data with nonzero `db2` is outside the supported equations below. + +The required checks are: + +```math +\begin{aligned} + &P^{\mathrm{rate}} > 0 \\ + &S^{\mathrm{pmech}} > 0 \\ + &R_{\mathrm{perm}} > 0 \\ + &R_{\mathrm{temp}} > 0 \\ + &T_r > 0 \\ + &T_f > 0,\quad T_g \ge 0 \\ + &T_w > 0 \\ + &T_n \ge 0,\quad T_{np} \ge 0,\quad T_{np}=0 \Rightarrow T_n=0 \\ + &V_{\mathrm{elm}} > 0 \\ + &G^{\min} \le G^{\max} \\ + &A_t > 0 \\ + &D_{\mathrm{turb}} \ge 0 \\ + &q_{\mathrm{NL}} \ge 0 \\ + &D_{\omega} \ge 0 \\ + &H_{\mathrm{dam}} > 0 \\ + &G_V^{(0)} < G_V^{(1)} < \cdots < G_V^{(5)} \\ + &0 \le P_{\mathrm{GV}}^{(0)} \le P_{\mathrm{GV}}^{(1)} \le \cdots \le P_{\mathrm{GV}}^{(5)} +\end{aligned} +``` + +Initialization also requires $N_{\mathrm{GV}}^{-1}$ to be single-valued at the solved initial operating point. + +### Model Derived Parameters + +```math +\begin{aligned} + k_n &= + \begin{cases} + T_n/T_{np} & T_{np} > 0 \\ + 0 & T_{np} = 0 + \end{cases} \\ + S_{\mathrm{hygov}}^{\mathrm{base}} + &= P^{\mathrm{rate}} +\end{aligned} +``` + +Here $S^{\mathrm{pmech}}$ is the `mva` parameter. The `pmech` signal uses +that base, while HYGOV turbine equations use the HYGOV component base: + +```math +\begin{aligned} + P^{\mathrm{hygov}} + &= P^{\mathrm{pmech}}\dfrac{S^{\mathrm{pmech}}}{S_{\mathrm{hygov}}^{\mathrm{base}}} \\ + P^{\mathrm{pmech}} + &= P^{\mathrm{hygov}}\dfrac{S_{\mathrm{hygov}}^{\mathrm{base}}}{S^{\mathrm{pmech}}} +\end{aligned} +``` + +The nonlinear gate-to-power function uses GridKit's smooth [Linear Segment](../../../../CommonMath.md#derived-functions) helper: + +```math +\begin{aligned} + N_{\mathrm{GV}}(x) + &= + P_{\mathrm{GV}}^{(0)} + + \sum_{k\in\{0,\ldots,4\}} + \text{linseg}\!\left( + x;\, + G_V^{(k)},\, + G_V^{(k+1)},\, + P_{\mathrm{GV}}^{(k+1)} - P_{\mathrm{GV}}^{(k)} + \right) +\end{aligned} +``` + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------|--------|-------------------------------------|------ +$x_n$ | [p.u.] | Speed lead-lag denominator state | Not circled in Fig. 1; realizes the `Tn`/`Tnp` block +$x_f$ | [p.u.] | Governor error filter output | State 1 in Fig. 1 +$c$ | [p.u.] | Desired-gate position | State 2 in Fig. 1 +$g$ | [p.u.] | Gate position | State 3 in Fig. 1; algebraic when $T_g = 0$ +$q$ | [p.u.] | Turbine flow | State 4 in Fig. 1 + +#### Algebraic + +Symbol | Units | Description | Note +----------------------------------|----------|-------------------------------------|------ +$\omega_{\mathrm{db}}$ | [p.u.] | Type 1 deadbanded speed deviation | Defined by CommonMath `deadband1` +$y_{\omega}$ | [p.u.] | Lead-lag-conditioned speed signal | Output of the speed lead-lag +$e_f$ | [p.u.] | Governor error into the filter | Reference path less conditioned speed and permanent-droop feedback +$f_c$ | [p.u./s] | Desired-gate derivative target | Before rate and position limits +$r_c$ | [p.u./s] | Rate-limited desired-gate derivative target | Limited by $\pm V_{\mathrm{elm}}$ +$P_{\mathrm{GV}}$ | [p.u.] | Nonlinear gate-to-power curve output | $N_{\mathrm{GV}}(g)$ +$P_{\mathrm{GV}}^{\mathrm{safe}}$ | [p.u.] | Safe gate-to-power value | Lower bounded by 0.01 +$H$ | [p.u.] | Turbine head | $(q/P_{\mathrm{GV}}^{\mathrm{safe}})^2$ +$P_{\mathrm{mech}}$ | [p.u.] | Mechanical power to generator | `mva` signal; converted at the HYGOV boundary + +### External Variables + +#### Differential +None. + +#### Algebraic + +Symbol | Units | Description | Note +--------------------------------|--------|--------------------------------|------ +$\omega$ | [p.u.] | Machine speed deviation | Read from a machine model +$P_{\mathrm{ref}}$ | [p.u.] | Active-power/load reference | HYGOV component base; external setpoint or constant parameter; source label: `Pref` +$P_{\mathrm{aux}}$ | [p.u.] | Auxiliary power input | HYGOV component base; optional, defaults to zero; source label: `Paux` + +## Model Equations + +The desired-gate transfer block before limits is: + +```math +\begin{aligned} + \dfrac{c(s)}{x_f(s)} + &= \dfrac{1+sT_r}{R_{\mathrm{temp}}T_rs} +\end{aligned} +``` + +The realization below forms the derivative target $f_c$, limits it to $r_c$, and applies anti-windup at the desired-gate limits. + +### Differential Equations + +```math +\begin{aligned} + 0 &= -T_{np}\dot x_n - x_n + \omega_{\mathrm{db}} \\ + 0 &= -T_f\dot x_f - x_f + e_f \\ + 0 &= -R_{\mathrm{temp}}T_r f_c + x_f + T_r\dot x_f \\ + 0 &= + -\dot c + + \text{antiwindup}\!\left( + c, + r_c, + G^{\min}, + G^{\max} + \right) \\ + 0 &= -T_g\dot g - g + c \\ + 0 &= -T_w\dot q + H_{\mathrm{dam}} - H +\end{aligned} +``` + +CommonMath defines the [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) target and smooth approximation. + +### Algebraic Equations + +```math +\begin{aligned} + 0 &= -\omega_{\mathrm{db}} + + \text{deadband1}\!\left(\omega, -D_{\omega}, D_{\omega}\right) \\ + 0 &= -y_{\omega} + + x_n + + k_n\left(\omega_{\mathrm{db}} - x_n\right) \\ + 0 &= -e_f + P_{\mathrm{ref}} + P_{\mathrm{aux}} - y_{\omega} - R_{\mathrm{perm}}c \\ + 0 &= -r_c + \text{clamp}\!\left(f_c, -V_{\mathrm{elm}}, V_{\mathrm{elm}}\right) \\[0.5ex] + 0 &= -P_{\mathrm{GV}} + N_{\mathrm{GV}}(g) \\ + 0 &= -P_{\mathrm{GV}}^{\mathrm{safe}} + \text{max}\!\left(P_{\mathrm{GV}}, 0.01\right) \\ + 0 &= -H + \left(\dfrac{q}{P_{\mathrm{GV}}^{\mathrm{safe}}}\right)^2 \\ + 0 &= -P_{\mathrm{mech}}^{\mathrm{pmech}} + + \left(A_tH(q - q_{\mathrm{NL}}) - D_{\mathrm{turb}}\omega g\right) + \dfrac{S_{\mathrm{hygov}}^{\mathrm{base}}}{S^{\mathrm{pmech}}} +\end{aligned} +``` + +The $P_{\mathrm{GV}}^{\mathrm{safe}}$ lower bound protects the head divider near a closed gate. A physically consistent operating point satisfies $P_{\mathrm{GV}}>0.01$ so the lower bound is inactive. + +CommonMath defines the helper targets and smooth approximations for [max, clamp, deadband1, and linseg](../../../../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. If optional signals are not connected, use: + +```math +\begin{aligned} + \omega_0 &= 0 \\ + P_{\mathrm{aux},0} &= 0 +\end{aligned} +``` + +The connected power-source model initializes mechanical power first; HYGOV reads +the `pmech` signal $P_{\mathrm{mech},0}^{\mathrm{pmech}}$ and converts it to +the component-base value $P_{\mathrm{mech},0}^{\mathrm{hygov}}$. At synchronous +speed the damping term vanishes: + +```math +\begin{aligned} + P_{\mathrm{mech},0}^{\mathrm{hygov}} + &= P_{\mathrm{mech},0}^{\mathrm{pmech}} + \dfrac{S^{\mathrm{pmech}}}{S_{\mathrm{hygov}}^{\mathrm{base}}} \\ + H_0 &= H_{\mathrm{dam}} \\ + q_0 &= q_{\mathrm{NL}} + \dfrac{P_{\mathrm{mech},0}^{\mathrm{hygov}}}{A_tH_0} \\ + P_{\mathrm{GV},0} &= \dfrac{q_0}{\sqrt{H_0}} \\ + P_{\mathrm{GV},0}^{\mathrm{safe}} &= \text{max}\!\left(P_{\mathrm{GV},0}, 0.01\right) \\ + g_0 &= N_{\mathrm{GV}}^{-1}\!\left(P_{\mathrm{GV},0}\right) \\ + c_0 &= g_0 +\end{aligned} +``` + +Then evaluate the governor chain. With $\dot c_0=0$ and $\dot x_{f,0}=0$, the temporary-droop block requires $x_{f,0}=0$: + +```math +\begin{aligned} + \omega_{\mathrm{db},0} &= \text{deadband1}\!\left(\omega_0, -D_{\omega}, D_{\omega}\right) \\ + x_{n,0} &= \omega_{\mathrm{db},0} \\ + y_{\omega,0} &= \omega_{\mathrm{db},0} \\ + x_{f,0} &= 0 \\ + e_{f,0} &= x_{f,0} \\ + f_{c,0} &= 0 \\ + r_{c,0} &= 0 \\ + P_{\mathrm{ref},0} &= e_{f,0} - P_{\mathrm{aux},0} + y_{\omega,0} + R_{\mathrm{perm}}c_0 +\end{aligned} +``` + +For an unsaturated zero-derivative start, require $G^{\min} \le c_0 \le G^{\max}$, $P_{\mathrm{GV},0}>0.01$, and a single-valued inverse $N_{\mathrm{GV}}^{-1}(P_{\mathrm{GV},0})$. Starts that bind a gate limit or cannot be served at the available head are outside these closed-form initialization formulas. + +## Model Outputs + +Output | Units | Description | Note +---------------|--------|-------------------------------------|------ +`pmech` | [p.u.] | Mechanical power output | `mva` signal; oriented as mechanical input to the machine +`leadlag` | [p.u.] | Speed lead-lag denominator state | Not circled in Fig. 1 +`filter` | [p.u.] | Governor error filter output | State 1 +`desiredgate` | [p.u.] | Desired-gate position | State 2 +`gate` | [p.u.] | Gate position | State 3 +`flow` | [p.u.] | Turbine flow | State 4 +`head` | [p.u.] | Turbine head | $(q/P_{\mathrm{GV}}^{\mathrm{safe}})^2$ +`pgv` | [p.u.] | Nonlinear gate-to-power curve output | diff --git a/GridKit/Model/PhasorDynamics/Governor/README.md b/GridKit/Model/PhasorDynamics/Governor/README.md index 8e6305aa8..0a107bcf7 100644 --- a/GridKit/Model/PhasorDynamics/Governor/README.md +++ b/GridKit/Model/PhasorDynamics/Governor/README.md @@ -9,3 +9,4 @@ A governor models the control system that regulates the output power of a machin There are a few standard Governor models - Turbine Governor (See [TGOV1](Tgov1/README.md)) +- Hydro Turbine Governor (See [HYGOV](HYGOV/README.md)) diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 32df36cdc..e51950e9e 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` `Tgov1 ` | the TGOV1 governor model | `pmech`, `speed` | `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` + `Hygov` | the HYGOV hydro turbine-governor model | `pmech`, `speed`\*, `pref`\*, `paux`\* | `Trate`, `mva`, `Rperm`, `Rtemp`, `Tr`, `Tf`, `Tg`, `Velm`, `Gmax`, `Gmin`, `Tw`, `At`, `Dturb`, `Qnl`, `Tn`, `Tnp`, `db1`, `db2`, `Hdam`, `Gv0`-`Gv5`, `Pgv0`-`Pgv5` | `pmech`, `leadlag`, `filter`, `desiredgate`, `gate`, `flow`, `head`, `pgv` `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` `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index d7b4cc9d7..26761302a 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -297,6 +297,38 @@ namespace GridKit addComponent(gov); } + // Add HYGOV governors + for (const auto& govdata : data.hygov) + { + auto* gov = new Hygov(govdata); + + if (govdata.ports.contains(HygovData::Ports::speed)) + { + IdxT speed = govdata.ports.at(HygovData::Ports::speed); + gov->getSignals().template attachSignalNode(getSignal(speed)); + } + + if (govdata.ports.contains(HygovData::Ports::pmech)) + { + IdxT pmech = govdata.ports.at(HygovData::Ports::pmech); + gov->getSignals().template assignSignalNode(getSignal(pmech)); + } + + if (govdata.ports.contains(HygovData::Ports::pref)) + { + IdxT pref = govdata.ports.at(HygovData::Ports::pref); + gov->getSignals().template attachSignalNode(getSignal(pref)); + } + + if (govdata.ports.contains(HygovData::Ports::paux)) + { + IdxT paux = govdata.ports.at(HygovData::Ports::paux); + gov->getSignals().template attachSignalNode(getSignal(paux)); + } + + addComponent(gov); + } + for (const auto& excitedata : data.exciter) { IdxT bus_index = 0; diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index d4deef66e..017e512ca 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -41,6 +42,7 @@ namespace GridKit using BusToSignalAdapterDataT = BusToSignalAdapterData; using BusFaultDataT = BusFaultData; using Tgov1DataT = Governor::Tgov1Data; + using HygovDataT = Governor::HygovData; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; using IeeestDataT = Stabilizer::IeeestData; @@ -99,6 +101,7 @@ namespace GridKit std::vector load; ///< Loads within the model std::vector loadzip; ///< Loads within the model std::vector gov; ///< Governors within the model + std::vector hygov; ///< HYGOV governors within the model std::vector exciter; ///< Exciters within the model std::vector sexspti; ///< SEXS-PTI exciters within the model std::vector stabilizer; ///< Stabilizers within the model diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index 00a6278c2..198a75dc6 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -142,6 +142,12 @@ namespace GridKit raw_component.get_to(gov); sm.gov.push_back(gov); } + else if (kind == "Hygov") + { + typename SystemModelData::HygovDataT gov; + raw_component.get_to(gov); + sm.hygov.push_back(gov); + } else if (kind == "Ieeet1") { typename SystemModelData::Ieeet1DataT exciter; diff --git a/docs/Figures/PhasorDynamics/HYGOV_diagram.png b/docs/Figures/PhasorDynamics/HYGOV_diagram.png new file mode 100755 index 0000000000000000000000000000000000000000..fc7d7916d2887f6758718be182a29573e7599993 GIT binary patch literal 72381 zcmbrmWmHws8ZNwPqyzy$I;2}Vq(Qor+DJ=xhe%4dq;v|1G@DLAy1N?zkr3(p=JuR3 zzA^64dl@LOS$nM+Z#?hwypwP>l~>ppBp46~1Y2HCS_1+>CWAl_Swxb%NMrrux){P*?vMXt5~y=>p>kMzH*$&qmWd-0b0|Lw((@)#pd zWNEC)N0N62(?t~G_<3qZ@YGM}Y2Y6xxTs`d(5g&^+5EgOZ&rSD1<*{j>Y@7O3!fFL7I6-fIX zZ-t%=ubNGS2DuXf13iz?kMZAXz$611Pk#9q!xp@pUou-OY!Blw5zr^@?5uaU*=}hg zBOuU7=i^QDnv3H_RpM|pyfY%pYYwK#SwhFRri3C_ZJT8L|Fg<_YRWFm;c&tf!I!zE z43^lP*pbTm$6IUx^gw_A9BTsG#jsVPvJQu4Os5v!w3hfVXz1%g!D24$!F2K` zGAMzI9#Zd=7(Qxeb&u3!=8`#K@~Fmiz*2(quaAZ=H)nkaWL%|glj3)&AQ%z|64c_k zxH-5+nLi5E7^&t6d~kF9C3D67Hy)6l1X79THq6=Iz^Kr1krA+?lR}U~sH>H~wR{dM z{*kIU+L!JVk*BZig^!09pPAgeEGI20!^48RR3&V3eJ)#Jpg&%v=i%D^`^8(QGjS|Z zL4thQAJhfj9Cda}v9yhdRE)Se5Gq7-JX{onKu0bkBLW0;2qJQrnmDo%s!<383=KjR zfhdhE{t*!u5dxuSGC5@8F>bovUp` zh?0-1P*b7rIIdN+GUDcqyYx}LojJ#fskkE)`*4Bwg4@IR=DtSY89_j~L@?*Kk?$O( zb(@X|V|-07B^7I`?%u!jh1R1LMIhlV?IAk*>4ox^-hWRXMswUyqW=03gY3=RyB0)xJ$k>_l+pICX0|N8zs0M`DXp*QwVP@RtVeGTcMj=Re zUZfV!aCAxQP40TiE{!^O=3p08GY(4;8Tej|%uGLa6n(`f#w1QcYtg|bK51==sunZu_=;{wCN#L^)Cd*WP$K7vGgeH=Qgs zpAho$fNEZY4ub1ppnc|IzH*K~L_MHm4c#?jFRU*H4;>rJzZhNe=2GCT(?e^+QC-`v zqY>fLQNw_jnH61{v>}@d4Qd}zd0HHwlFG-}#;05{(}^aE4ivjsP%R0%KW-l3gAyzazTrk{7=HO|O3H>5{aF<3rLeB_VP! zZa~x{Q_w>LyDkG<-$pHqdCN6*=gVAsp z{~HZY7b~-oz2776=mwMvF)uNmN&2;w@Y(gR4SeXOB@&PUHD6nn8Frlz2Vaq>d>I6E z3=GIeadZqkTVhO0)vo|nFbuozhN)9wG(ue@5W4g zK#R8%h{Uv$JV>5?PIKZW$Z*1G)}z$&V$**v#3o`lFQgu@>cEms$T#3Tas-sy-&WAx4 z^rk}`kyCM^#Jft8Y!3FW=^Bdexh-+B4`(dz$T51&Dn!Ed-XnG1Xsy-o&zkmr)_MDI zYNw{Z`b4uEGG4@*gYj;W-hx61{ zlyo(p^=~bW^B>V*o&J_dAoR6_aXDn4J=Cmb#_Ar8bOm}Cb``w$p)JNGo_w_e;{3bN zXE@SZPb8@Gnzlq#DEoWNOVdM%`f<=1gNlB1E)~z;rrvr5cjW-gdfFtDCR=A@rH5M+yaC4sdtAyaJ%Sx zWUcg>3mIIxV_evGq%{ip{Ov->5UnQl-4X8dlDlD`Rr zxJ;X)cLEd=|IjBG&)|RvX%R7HHezDF8hA|HW`jaG*2}cEy zUaaXj;dYDT<#9pQ$5V1tH4MA!F~!k`IFZnwFsr*Wqxn^ceauJME{?K{tTPUG z8Y*-MGK#eN$C4f_ZUS^QWGXdkNo1p^2rmB5Flrz*{pjcc=vZpvoTR`-BnS|3K_v5+ z{HV@}B+DH|&yAdZ1SJ9++_yM~O6IJxwchid2jcIxYE7Jv#B-Y>yg#6|4>8!E&#E_0W7(@-6JP1&+U>D_}cGZq)=R(M3um@ zg2KXzMd9|eFO8Qsk5$u;m0-r5&SB0tWgdOzV)xwLvMslyGd z|DCg!8-_!b_nuJbVvKTge6}7HrHXXf*dX^8ljV9SqiH@D2x?#qqF_z*hfFTr1+f|# z$jgm#qoPBoFW!|&ej3o}MF^BxW$$HCRW9SDNQQJR&Pnu~a3r(z(iY=s+XL*@SbAzej_B4VWtVg5 zQ)`?H&0{cMf$Bfnp=aa;L&%F!!A#udLf8abI*3FAbC#^v?u{*^@WnFwE^tcbT{F2u zMML9^0>}Ep&zbq2s8_{9{Pm97@bbYw?$ysRkBQn5AUenQ+H+M_^}p4nr=WD$1H9iB8(VxezoN*+!Et_AHlYFaGW<>;o9yZ)0)G;F?qNOk9>>>Ur zL(WurvT|XOoxr-_DPYKZ^ol3!o3b94)pv5W9$CGIBnv^1FeBEb(2T*Nf)RU_`GU-# z4KGu+P|Dlt#2|1eMWV!f%xEU~T;C=e+qm!LBb{nekVHa(UF`E&eHUnEMJp+n0$2h5 zEO;#bcJWjA+Xbj}-cSP54{Pc8RcM2tScV4f%y#WA_$dUZZpJH>XMKYJ3j&Xdz)}DC zMl9Q6JZx{n6;n^Jks%(;*1TPq^2qnwk372bc`UHnawJ0N`)sdy!@|NyBym5KO$p;Q z4*xHGofAvgHWFJjkE!e~eh6!sH$8u9EBPn5=ENMng(ohUm9AZg(m!Nr#OfiKs}C5Z zDdkdcY@a9F>rna0<^LbarKP_Jx#u*VAHdvskh_Zsg4V3I&W{#6T(@tM}*g7S!m0w~HXpG{URQEut2- zT{XiFF%bz zUtaD~kcLky&KeJ#-i*JjBXzf3h<@FO0)KVfV>Jx^oAnz* z_vfrVQO|{T2^l<%0#2yTG6w-GfBT(()hDPZ+n1zNtL%G{wgJ--D5g)mhlXis+VkHq$7g0@ z2X*>yZM`GLLcPMRphJ0FVAS=ia6c~l{lcep$EN;C&agg@Z~DzoXs{|#9&r|pw zr1%K#8nJ)~Ak25rqd9wCZ-3Q@a$+pKWpRGA8tb+{v))7MkmvzS?+Qz$#4?zfnK|zB zXC?>w$9oGf7SzWStPWqyO@vz|OQ$6Zt+T99qcb%SAhRqEHF|t53m+bON|q+U%n&jL z*gCfMF7$c@^$r3-qVQP^rSwHF(~3-D(WVg+g06s)cnhD{@!u1Gg@6H~dl5D2dU*-* zIjju%OJUlC=SZ0%;b=I?F9v=BgPh5W*Qk_$KEa2SH`V!-3E{xauFh$*Cs(_|!1Lyy zep?^xt!NNJsFoatz_H|9&86u3l%F_4JbY_w?C^`rzkOqkdd7{an8N;}SUrDh`QP=t zz|J)*%l2}{Epuk~x1~VL_}=8U?|64MYr54Jf5HxeQucRXF)$J=0I2F&SXpb!Rq%+3 z8UJkw6S%W_>GvE^)6{_F>onqYbi)8!h52oL45O+}{BPXK&a9yy`OEj`wRNwgrHxK^ z$K6d!u6R9^Lqo>)&*W58LbOVQp=HS)`0{$wSqxwW=ASuh^oK$qr`|7Lym*l_ZN04@oj1E5pOFz&Tgy}M<;!;?&M6FUOP)07 zs%oKzM9G^4*j2=IskTeuW14pSWs~3+F|_t<0G6^!zRnswhv(allhyYRf>>UUwdJp2Dur6 z%?i{gIi3Gd%D&Q+LNBM}N|zkPQqDZtK5$gj!Hk^nS~ueBAPQ(c!az@9oBNXK5%Fk~EdrkpQwp zNB>gffi;W-JkyuZ^LWdti#WgvQ=7!(IkPE?pNVN-H)V z5^TuS03Yovn39sJ>MJOeB+q1n>+f)^B@jxOdFv%F$GoCP%89Gw*v>yC1v@lK}+ ztGS~k?AHTC%08)rryaj6HZJ_aBhba!Jp}dN@~>alltNS zh4sdc-#W@1OI5>Z^|>_ymVjDmj}cO1gaq0N1M~C`tTZkU6TgIjJ(QOt^UM`qnk}v$ep5OcJGbT3t&hR&W)4w3qm?A+XzfDvr3X;ut zx>Bi@Zpvpppj9x-^Fdb+vf?56HafAn>M->(IFEgz&Fcqm*YR2b0Y!6I8lH_ntz=et z=-Q)0-{E0ATgz|~9F4C=iX*v6TDFeC!aG5w)7HNi>mOe%xy0(9sHwjIN-JUwRmt-@tnfN47J^Uc*)S;M8y>ewjCaY`m8 z>s{A--OY}48-d+3ZxE|19)^k?vJ=mHnN(Zm)`n1DfxB(~L3>3@|5dooZkYz;pT~eq zg9*W>W3Tbc2JpeHL}gjYAk4HH#1<4hFDNLWVq$tEl!>EI_yq2KopB*=85_S_^5nI& zEZP4>_OK9QOd9Z3#n68yb(lN5p9_-UxaTank2*nmt>6Bc_pYmulGnx$dK6Nl*My=k zUzSb3{dN<19h1=C-{0VUTnYch{44V|F3M-mgeqF7E-Xl*QNiy-Yo;9aZ{Ar$Ny;fi zhS)k$b=q*)i(ZO@{7=yVHlO*P*w%j=CxAb5%*?_yUI*fkME7h*%(Y zApV@mJ%DO#^xrL_J}F2*@!|`rAByxqvyecuDsQ>HD~YUOMY9M%r_nYuGi&6DHUY*8 zk8WGTD!{F2sW@?)J&t-J2wAq4x1`=SHY8Ykz-_gn^*B{Aprkm}v%BY_gVOIde0wlu zusJ+qV!!6s5+y>;d)EwDs@)qUkvfNaebMgRQq;#?jWrfVoWwdGKlgQkdZfMADCj*M z1A|n>tX`4ga|VXNP#yV6rd#UtT=tV3DQD;Ev#n(#&IU49$HeW9->A-^!BBG9RbKU# zw~;I!_1ma8&tj@iFz4F+#M}>;06eD59sT;SVSjdEXk2Z!9Z@`o%<9S|;5Hpcs`geR z5TIZqkZ1rH35_wcpk?Nh00s~>>~!#pntyGoAK(3MWDeCRmP~M>IGsa!iv*RduYzlC z%|Jso$c`bZODc%b zW=}M*M7+uC1J|69Y8t9>vpi96%mR5l`^UoXCT`0gQnHG8^*M!0 zXQ>{GPyOH<%q_o2>C>E8{i}Jl)N0SVM*}M^Mn^$GfoJ;1)P__MX9i32=Hk*nx|?~F z(fB+PA|oQy3=BqS^0z*H{D_?L;webe{#K7;Bz)~-26@RKYp%XoXAdJI@574s;G5NL z6lAI~jB1r2YP5YY0*XgV7!*V(0G<%tlv3ipysk~gG%qm^qK&Dh_(fwTSo_5Y+bgDs zJE|bEABS8BlB{Wc8!%G6R9P>ls#?!J3tgR=i!8c+L1tTaf6tQgWluc04M!0}Sp5UY z+&DxwFP6M_korc2W2RUA=#HG_21epz^ejd);-jf(JODDh#R;8RvsOjHyW7BMH}(W6J$7;kW=>BMHe zIIEbp7%q!cv+HgaW_Qo7=?xp2(?PfgdH>U=PrqX(fsZpWF^yTnX85cY^rU_;FjjT_ zV0r>^WIk#b+}l5i6fVaQFDH)w=FDWe{;RW#6J${053|xKkt+(lQdZeK>B9|<>w&fe zc#TL%I9u8Px>Ol;A|j)qanG%pPw-&q9CnnxVJ;%-Dzg&;HGr(UJ0Hlm8Rk9$-W3Ho z(NySiN=niIXuBpSu|-5g#DRH$nhz{H^d_8yR#YN9Q|}>;zgf9xB(v;+gNdnX|1U3d z5(dar>O=>}YS!+uvZ&QHHD!ueEJ+FgQUKTrjTKZ-)X(gp@F2u!@ZlE(XsM}LAQKTM zP=9A6v*sqPFCL@~RJKkbvjt38wN=J^_*Gh7GVFgH#9nMa(B`y}tk-{geegf|DH!UKO z1D|hOHrZAG&ey5OrWV_~+Zh`hyBp%izX!#QuS|prW+2tki4EB1XV0FgmAsjm^3hh| z0~j|hPVFH>->S)9Hos~mSi9s6OvU}o`z1hH#sCNpn$Pf0cO=d3*<_$lf@%S5i+!0N zH46T(Ch@pKbA?e0QqSYwixE<}&kU8wqhw}9muVmy+pgX=F9zYjL8R$Er8`j4?Pp4?oN7s@d-<}rdoN+^?l zQ*p|gcLP81Y=tRNWl2AC@L>TnHurVWi2KcpBNP6pH{6m!{yC^wl;w5x>yhHS&%f_XjpP6Q^MF#?S>34lcM$2oV!+2vIge*`{W z?Er~RK}O)q_dSRb6>(?fje)gwX(BeCl(Fii$;S zL=*=6y_)pZR+L5wywKS56}Yr7Y2Uz2 zyAWnBa7kh{fWqp*k&&r$6M1^-3P5H33R^w4Q1P^N3cF>oTf$)u$K^40sMRq@KEy7w z()(Jo{mrh_UelV|(?U*1)1@iPv-b)`lhe=`jfz=e+X|Z`5~hy0w7ksNhp=>ke56+{8Z^R$k z>-=82m+_xo^p|r836Z22c$est;x#_580f3u?CQMPa~0Z8d<^i>MKE1qUR>R0?1^9a z(eqmYC?%;&*!**=X7tX_K^%GXNI`_lDyj=pLECn<8sRm=N(4p9X*apYlJV?pvxkFg zoo?0lY+xk&r;Jt!TH$f^vxH=7P_2|W_VfKT7Y;K;5EA>Fm+0x}P+oGVZ=4Yo@p)6m zA3k>gcdK9JnlCa{s{Wy%B(t>b=1J^>X>FG++N%Pg#~S8c1!xb{lpA%V z2+~``@Ccu}zDBwgJaYS)hkx<>Be?Mo=$K+!WSK?@nf-)_oz*vyH7XQ)z*M~S#*@77 zGV1H;!|P8yUn6Esy)h3ue9E4BL;vNa04km_HIC|1D7H%$p#*GRB|q*(&M)R#mS~T{ zIR{gwB!zdouBk6{qk)-EZDR_DyD|vw%t_cd&Og1tvt|Vp@9*EXZ@G_|S;77e`QGD! zyVGlZxF}q9QcCl!GCE8r7RwXkD|;YhPWq~4j>LLnYSTYJW6jh8eBwJ6+L6V3L}q2Q zV{M^9YielOOu2qjLmaRlkQj6@d*6Ubq5d z=!#vn=DKP0Gr38|?(?W#xr6)juVxM|w~aEL<8@LHYEyakB(6B4h`GL3{Nm>+xR*4Q z3wEKa>y>$hl}s2mhUM`tNr8^uX3BGGoeBOHFkrvdd7KHE6=3r@roAejCJUReKy1LP zb;D*2Y-VMra6yqHR(^@F>U+a>f{pfVlC>vxH&!^|daY@2!Spul9)6$8!$H*?jgsT; z+FtVeVJtEMT>cf$e4y(tDd86CeV?)h{0kt>5Yzb`de%CFwwDe1r)M*=vtxI6Eq|}J zQ?|apNcP}o+2a+vJ??U6SO>h+_ucXQgC3!LtN9 zP@2FrTJ{yp&TG57rT1IBUV%JPru&FFX`ty-F!$cxo}8SVqo;%kNG+=zFJF@X`QxVj z>zuDtyQ1ggBYXhV*pI%>Zv2ScRyH_lyE#nYpO!D0T%Whd{!wc+`Fo*|i8mVNabzHJ zIwmFHv&A@)`4S1zV13Lf-Ag|Om38UxM*s_sqn%4E;7HXI15?esK#7T|0z8ARuI`L@ z5YNHj_-bZGV(`|0&SD8ZKuBuLNqqwYU)@^bVy@uV@L5#Idk=%V5G7uD040D zi{CgLipDJ=Ia+FsxjqNPGh>_D`Xn`~3U&6&a6T=+xVW!LzL(*ujl-i0o8PEA3=dG6 zCx8`)QfqUpw1$^fTy6CygoTH91|Xq1&YPw5H>YcUKnOiOJ@pn-BhSPt0>seom~{S6 zKlx^hm>AfwLCyx$+xo;`j6i?j7)1=Cs6n0fNsnMq{hjbq*YWk=?->{vHoL>|J%MQP zr;j$auUStCz%GCSf8D<1qg&&&nPE*#O6pHxhV3if zw(%-Ug=OVOP6x2rjbuIhc>a2UXxX5tsVO1D=6q1|e~M9_aE*+Lih@cE^-wSX1b+Od zm(uj_&tbsZsv#riH8c?CG3I?}u7h)QWl4;-GOu2l0@_aH#xAJhQ$vFQcuw2MYBob- z^+N;88k~9PPj7)r#TvR<{T5{sbo)gB$V(gw5ym5DMs{`_V5zdMuCAHwdNdw@lJXeb zfmN7vZV#Vy{%KM`qQzfsR{tok- z?vCZ2P_s$qkBb*gN?38L>X$PigwAv%AfWYEwBJWWlkyLwXZ{+?MUsso4!FO&8EkF| zq(Yw&-X38myE`D(Gt|<29g|>Q^9=3+hJJT+M_(s%B!W1trx=T1Ng-Iz0^1+%Ni{V! zPuN|5O+6Y+W^?fzI_A+nr}RJ}gVuccoltl3_ivTMswjY+1Ox;*O-&>aTwGlLa{Xo+ zQBex;&QxxTm@|P{n2k120jRjSC)bZs`t%Dn<@I-N?}z@3XZO95C-PV^=#mf`BW9iM zNk~Z0)BhP*k>KkQ9v+TUI$&XDCUGDMg=v)3He;4d?q#TaE;&L{D0!2YlLL{Fk(sNt zN@!{RA>Ww=^`X^FmBK}?aU zlxU+tzFnhAr5ZqLQ+;U9%cDt5OkC>-Ksuo|2nawsnl4l%T5=v_OBMYfbbnUix88|H z*$v{DX76&NQur}V=o zA&D3pQ|ovVl~AotCI~dS!}bSXFYfqPwLP1Uir&;pmaC{bnoob#bobyZco4RBDQpLh z7f6-!;pSYWS4&gVc&b?aayLKL@oGCom<$O5%EKPoSdl=uz&SoX4!gE{I4T~Ve*mJw zJ3E)88YaTR!a@O$!?&4=+45Ey8hA>6@)u{%Us-+5V_YDFx}ZtD)Z8yj6wFi!otk=+ zQcd=}?r&;R@^(83(AF?2Gi+gYi%)^dnf%zwy}ZRDR`i?d)lR%sQ8L#}%hKgySvk6KILK6;KG zXvFq)&#G>!3Hn=dCM&g(%<|mr<{!p;FaS!{dy=B;vhlYx(W8tS|R7;tla;6X=h+YZ& zvxB-F#``m6Ixy>`Wt~?bewdh;IPMju)TJe?e{p>3*w@!r>8n30Qy~nHi|g(f!o&4! zI|+fKKkF_K7w6X+2jk=8896v!IXZG8A|j5}+38*33J;%UO<8PohhN^`96fV%lmscS zs+Xda@GX2Qh)75W3%@qm+W#igi}lP95&uO4RU8Nvyg>FkSjZ%Fu|@auJKBN|#?*AY z#zGa9EDh)?iyu_#3xI*R9VkmN(TMy-yZKumDBky zkGmd|x{26gSX2zVKm2(!F6So_A~;Ewy*%QZ;_R{dfj^m5m&SqK`h~?!&cUDVY%VZA zFTDVqhn*%684%E3#qafGgpshqu-`Ci22ofj;vy?$i!OWjo&brh!cDuS zQ4m9*PkYnxRR_H4fBCXiJs_(GmsPXK4j93Ho5Ad-V(cIOQeBOk#QP3u`vQAkL5ZXr=0T#apg>Xi^Xno z1)D0bzvf;w zz|hNb*at{Y>i_W(i$d7qrPDOr-1tD1jBYbM8cs|$i3J!J=F7L3z^5(G;1~UJJ)ZrS z-DJKpeZp{9GQ$A__WZ0qpv}{H=?p-oE7PJC7G8wf(DL$5H;Hx7U&^wTR`u$fw18s% z^k|h*LPFw%y%c0X=*rZ$G@A#ja?=#ReD->`e2R!x$eslchKWx4%;(37LTyOmGv;mW zlk-Zans{gbbT8Jfl!}~?coOh4F>lg^ZY#Xrd<>dy0u6Iz`X;IsV%GWy7+a9@A%$mX=#-q7QqP!!+HW$6q4tl z5RhQ#E5A4RMS~JjU=&)*A;&;4RSn$_&*T8LxlWU7iQ6x2OH0PUpdbVYsB@#gnr^ZW zw>|(e@Ht@V(=)Aq?RZK`N|s${VSM5{g?S1mDFg$UIHP+cTc$PZEf%&^*ci*O5`utc z?Q}Te7?=IKbo!T2ij*}Im^YjASbtxiqfG#PP^a2{OiQOxvD=>_jgs4quD696)6Xp9 z7};2&8=N&5*@E^IDuGD}x9S*x=(l{D2|99wR4daWwMy;cKeV#fDW>v;9Ne?r#n^F^ z+|o7t9E4AzLLK_007wj*X*T5&Zf?A7_t!RCKYz}=&p3PTu;Pb>DwPmxTz6y#Uzkzb ztJyZCk~bTR4(&hwj*7(s3R6dnsL*A!wzd|!!zF#VHQbxz1H_?8SkoLhKi1^a3w*K38{BFf!t!WXmT$uJv|Zu z2~}rtOOa8dOa9qwAf14OjuZ4UMrB5C^Im0&06=a~;2^WQ9O}`klJOTnu_v6AOm13% z{F+2c@ z3&xIGYz)-K+-iaiMR|Fk+q-y6(wM1Uet++Kye6^aHmiST4~JEuP0|?7f|yT0%?J8$ zQYJ$ROC(}eZc+?=8)sa)dwUl@C}XDpG7@Aozz^-ZHbIZLSleHD`A5sqEJwSWB_P&G z9|E=?u6To{4?uo}KJ%Cm-PxV{8T%%A&qKe-RRWBZ+&Ak<{QJfu=GbDqA#?7v9B&&b zd5_gX;kdeS>r54kR;g^D`zB4Yw0w4bHOuExP;L}xP*G8{nhtUuO1AH{RjZz#$w6Qhm8e_5vleS9);do{(TuhD;+51WVxLW%xY?t={us&nnzA)cC1QmF zDt13uUs5vA4Fkm%<2k7ICG1KV9o{3XAHi>okAJ`ivOe5q4o)yl!Y8Vgk_m>?XR0Cs zU5mQ456^*cF<)556ihU24THsm)vi8c9W?w-59oY00IjcY;tfv~+EV~PT!8iV^ejxf zSRR1N@!;G8!0!uKjd%gbAI{D%fW>vy>g{<9B)B=t0y(Ah9&pN_(E}m5-E8mjYxPt^ z$$nrQv>u2;(Cqn>QFlQmT&VWXToIPUu)bvyzob`NsLI2RaB&wXTW3<5_nn{P z4~|tBiZ$2GsI;{Y@(V5w1vjb8ow3h2IG*a9*m)~Ts9ObnGYIPDgVK0a^sdR*oNjmeazr^A$1%5|3WXrlB7%kbHhDF{7EGB* zp0KZAR}ssD%>sr7q~0;TUy`_jCZtAk%F1pHstS|uw|e<#XgomEDfG?j7i0h!VI}kU zwHC4mX_R~%sSXT0g%Nnd)tPi(g!IUFhGy>(*LUCLF}Ws7hEBYHj^ebn*Qt^1iZGmK zF~Yx9 z0dGM1$GDXd#kzfH2Z{{ntKs{Fn9!2#+5!R@Xos4q^{@!A+w=UjD*z?3602<$PxQ-* zKPo+{$)P^=ZaupKQoJ>JP=O2Wt87|n>BFhxQV|;d>y1c`Ng=1DVdE(M; zhwCeF0}Bg1q{m=A097>6pU&~=Fi(C^UtS{RT8?mhpLRJulx~ummbOdCVJev*RSw2V zo*Eaoc3dD&VnLeWuRJaV08q%x8mQEJe@&N;xBL0|fW|PmtGsZ>o2^@1f`Ym-QCOM4kW+~g|;+ma_USB+dWYH+898|K3GliX&HI(K{3sGYu7>5 zqNHlQIQBMZr25U*gOkA%(}IHXHO&u+3WaWFn0q6>djh%B*80~XDYKCHNoyFO zW+c76o7|YQ+4fxY29o0wVoD|~`0>_m4X!fhcI=eDId<6y>!#k52v-`2XX$PIN(BNn0O#1Cacv=nd-^ z<#h8|cEaa{Jyru;T1}TJ5r1bGUcOw0R@QTUel&6aoG`*}$Reell7ScAbdlM$R>xm` zcYDGQ`VWAz1VUr=+1~|xsCD-5$)d*u7%1nz=V?L58^0-Nng|O)ChZdg`K|a3Ga$s9 zc~thwU2UemD3dPvJ=}S&1)wTDdo%Ti#-Me~eXHw9RN+A)S2=0OwH$mUXVpA8iTc`#s5Y<+S!ny@ql5SqCee_?v! zgWh=%)ZDT7YDI+<^|1>JdqNUBjc)t|bTkWK#5|U`=0m9gi}uafU?luTUK~IP*pdPV z9=uO?1Gr3k*7l0Bdn%2)wxs$Zc}_Pz`GmUOI&YhNt5QVWXFCght}v*;?}kDXo(9MW z8U_Y5Z$tF)r$V%5*# zt_iRurae)eb5*8X#@~@e0MY3Qs=4E}pvN0<^kig|v%J~gEa&m1|3U7V-q~<|eI`;x z^|zeFtgJ&$VTrGQsO}7}wLv3CO6Sy`YhQoA)mg8?iO8!PMqr~ro*3>~48%qaLXel2 zzg4qC)ga-tO?G`<@&~ep5wL4LWiHbuOCi>x;N=T?loux_b zFE=q(*$cfQSllnZDUTtI(rFMl;ZmYK+|R4=LGA_V-?)moW{A5z1{+ z4!SWTlQN@y?Q)>ko(f4_F|PV?;ayzzs5Ja|ZCtOHX=MP^0QAOjE7=Mt*8HrFXnUm) z#Nez6@WsWm6AYlof{D&3k)o=E7ncAvJ=UZ$AR~Sv=(DkNeEx6LM>S+Vm=`wIjF~FE zn8GI{6K|*v6Qs95e`{PX=bZQ!iN!qj6LUe_D~_N}h&sXcFz{{P00guXbX3x?x7yn7 zK_7K-r0>J$IYGbZD7b$;7*Me>P+mcou-&qeyw&i?2z6`of|R(OFqG3dJG+{@;eh_pEbkhqFs$1+&0_9@Em(_biB7j!{KX zyj-yV#tmWzptk^jk>n9S1OjH2%`E$pQ!lyuL7UEvI4&}9%XjAH(8aG`zpAOJ8S#&N z1`McI!mM7aH!r=)~>0cZyJGTy^06 z3hzGq3|elB5oy5MY!4eaoc)?k<#)(xc?McC7+E%#6#GA1EL&XgBH<`cF7?i(rICP2 zBoAcD-MK0=72nn?fxSpjVFwtg9*Et4r8}1g+(m3DLTG!vZAVg2Xn*>02tsIbm%kLt zR=7XX^v@XoAlf`elbLV2C+5!ICmGnX))?2uplte01n(~nllT+B!H>h8e$(IQ2aA7p zvi&vyqjI6+-^SIlEwE>UfMx;jAXJ_b&st%4=4dHej!fozb-Ee~bex`MCx%=&gyN;QrzU7DFD_&Q7Z#nwY#yv~0S6*wWYO$|W?x4*j{_Pl7{8NC@T$o*}uibq23y-Z1NSA&7fRzKOpttPH-&0X=pW&i%5}bCQ2g3W6bB3(< zcI@z34}BKI*2bX-@o;0=l0s07GaebTO655dd@JZl+|=8xx9E3e%msQ|d%-ASb&Ft3 z|2CrAU+VAi8_<{nKXGwy_GjYpc!Ti1tyc-?~sD(WaYr%1Yu-!bFzdJfeSWlQz9idu-5?{w12ODn6>^@E6 z{NlfDHAC7L^2CE*7#1e-Za9rMEIRr}Upy_m<-Y8bNVCKQm`VN86};)UmqU!TuI{1f z)~4C%6BrIG5VuYO^90_m58Ac8#Ldj;fSwqcZr}9ipZDYA!4zVuNuO&p+a*iZ?+VNp zM{c=M^a=kv2NVoy8_*i=GnB?_QwdN7yfs`GoM()iGZ$p7o^p=s$rAGV?rZD%w#?`w z4`EOz`VqX_9h`ba13lyXOO8P;k)XXF%HxM)-?%B7F#v@3=Ixjpi0<)PvM^-_+$4v+ z;fr6Sl2yqb$`maIzGq6_h~FY_U^QpyHseZ1OZidIqoAaO3;|$rusN*n6`7zja{a2W zaF+KkCU|x{v3K?#G|_NbuTCBAb=MsNe(KnOD5J0OMbu(gD=xV;G5NdCA*>NRXQ_$# z$6m`Oe2VU;*xI16`b5Os!F{9^@ z8{(xyS=aTg}R$fkT>#)-Zna zeBh5rULO;eLEq!&ZD8hQJ>brc_qSef<3VU(WzHL4B>KD35u#~0Ywa)u@*po+pEU^8 zCtqZnx0ChaQc}<$sX}hdUsN)gDJQ<39V}|ENWc3*f+O8#Gf^Oqj)~dDknKZrl;>R; z6664nN`z0JhJ@oY6+Q7#^~E zz%~I1&nq}M*!cGH&o1cDfR#C{w)uJu;UOb1Cq*EjQT8$sZn=|0%dWs-ju(EsD(Fsr z4_t4YbrF;RF&P;d(1WanH_S=5y>!;R((2Obe0J~yILF71AA^Ib0k-Yz>;#9uB>^)H z$oL#JM(_!*YC1Z^kx@{n1q6se5=6(wmH}mI$YaVFYN-(ThMYJ^NlQ;Lb9DtQ4X6W5 z;=$UN5=}1r6sU@SnuXDUjI2}r?oplXJk1e}HLq5M0k#D0+XVs)nKHgbHXyj)U7scR zVLBs#h8%ouvNR%Yvc6e32mu-yI_QjdK!t$Pk(!rRb46(jzreo>uEXk!9sp}I2PZ{Z z5Rq`&fFJ#wUO)tWxg^jstrZ#3kFS~g>_FK#X@taz5ZZJN>t1%;$Zy-M+9L0utf|EjIFM) zmbAvp#8F!lv#Q%|@Fy^?=m0>n19jW7mZ317SLnoq#y5l31$$ax^VU`^?#CbC5fFw8 z?+pS03cLYyc>2VRTT1KdJ{!pvi*?xInJcVL+b) z`t3UCuIRFlhuj)9A>gM(?EzvpYJoVOJSo<4;OpQ|4($PN7N z+(2Q!5zP;r44YITuu&sq23`DbjQB+-YeV%dITW8@qnL$K9IN0{dBwY zH4NTozk=EIo2aO$Am%mmX6k!B1v?M+2L1p;pFb!oUOVJNSuaUJmo(SO$%zqP zI$W3j28Ce&$NL^n`^piVZqr{2N zxVdhYul^lc{m&VGlk?`k+g;&(s*eIY(%(t5`6Hk`#+SjHo}Ol8WtFzKX9vvyVL%h^ zu;kgez4{B!AqF4LG~pr(=xK$wPD7}G7L+OE1}Br?vig2`GY-%dU#Y3#UZ3rU{(o$} z2RN2*{60)FvXv-=XlM{+g^a9r84=2CT2aZ~86l&RCPkvGY$7u|iAa%6WUsP$&&&7s z`@jG9I9^A`caZ0K?)$#3&p1Eld0xA;-gNqK8Wuq{CwhDM+^vz8jTns)eG_C3%X4+^ zrih1|WNrzyWLX6QBsZQN?{XL~rWJFXym;cdG#mD`GDZeiF1{eXZ#aGEFMFi>vg3(J zxy=#hZ>+-*FDx|F11+ub=f`4%oe~`#&C1TcHvMN1s^Kn#(Nf3;{n2bwLA0lO?b%ezas!`^Y6L@@lWplI`S?dLF2*&qIL3*zFNCh+GSQWNKH_C_r}rO zeqQafXSZV2k>v@*2ZPP83A+1$fXARo#*|<5S?$%W$0E`RW`gN9W%4q!oQJM zXGDLi7MP$I7yx?l`|}gQF>`SxR&bXvGreZ@?LX-tKpLHdOzyFD*9(^vy{Bf@ zy!^u~Wb6L@4axT^pKYG{_tW(dSvP-2<#>M1?pfb$Iy}GakY9U{Ah&zCl=E}e^}Ze= zcLCC!Ugc|Sm6d%K{>Rn-|AK!EW35`V&gDCZLEkDkme!JEEBHv{LICzyW!|9nxpV$C z`wDp&IQE9x=UzB>Zc}*(UkKoKz1I1EBZ2>=!juPJ{4(^DBo-dHY)JX{>vdY%`kkkq z7iC$szQ6TH(*?``lpc)WeuGLzB!H-FKu73$myn!hXah9+#FPEh5azI`1@qj+mc)-k zLs^MD7@-Q|;FEzg+!W~LUZSETfC08D;1aY3wmtWxy0coXOFz!EGhiO|Uy{=sYrs2m zV6%VU#g&_zI}2oEd^B&Y8@*9E!ji_>v)fm0tGJ!6Hhy0VH9K^pMLWJhi!Rcy!h89n zEKAB~hhNv;Q&Jc@n!h8;#d8;C^iD^~Ghx*6`SS@CRn=EUg>E8N=xP+vAw#;{*r`ap zr0~&@MTs@h0Fb@$6H{{oV}rPFReZ+`gDL&g?(4V;2@CrO1yt75VEkHYaQ)+@dqIZC zIOvj*LW_Gs6oQE!cxARi!KN==FKnNdGDEq@nWDe^1t*+=FOh{r)5>%vp}lzX$EKR|$_AE16|@hqA27!!HyeUIV^ z1nzJeQ{!!p6OX7(OIvt0hr0nx8-QC^%bZf9y zz1u%J*F-%!j|m0OWyz`t*8aJ$F_WTKl-O8bA`b=PIdx%Y7gg}N2Pv#9Zote51)a)? z6N~`9eehi#r~&x|1O@RJ@lV@YiRPEM!RU|85}vBD^700|EzQ3|GCm{3?jIbyO3>gm za$(bYbD|g%dHe|~`6Jk)L z6jR@L@T~WFtq$!5%)uNt%z1pj+1GjmB(67aXe1;gkSV>48L=eEwF%x-tu17P`U@&#EebA-9A+qSI*-APt{wv)>8 zXBstL|J^%!EL$)pA{4fS)!nC8!{`ahV?_+wWo7f)Sx^SNNLMVWIUsT>me_ne*b#J- z>oLE;N;Y~Z2u4}6+in{qIjPzG%)PlYLthwOZ7QsHvDWoev(mI8Vf_@ ze2zz{A4?Au%w;SJi-bKV#7+G2-9!*R_`6y zh+VxJFVPd@zOw9qP&)qkKwT&3i^`Ukd)F!{Dk)KtJWQG!qTxkq%vszqC;ryW8G%2#@T2NXf%qDsZzksiVod(4?W|J1*tz3!LSc;@1q@g{0 zhwt_74fK?{=)y_G%nZ_A7y0meQu4iK+MQXrB`p<^Jy4U00*fL!3hTZc*Q(DFKXEutkJ+@4m0_yz)xrkQC5k%)l*V z*(7S=_Mu9=pgrTr4;;#Q$kT z%&v`K)8c)-{$qL#&ksgJS%_^42Awh3o0Hb&YyLd1d~_`K-Pw!JTxe4s3fG*Al(t<1 z@X5RW>L)Z;p)dq^P<02h!ox{AF(CYq8U|q8SCFw#wO2*?{u07EdiClR1zA)7hf9Sw z;RP)LX$;ASTgv(l7A#^y?`fspzm^Je;NYTsmRLX(*tb=roa+`bVyRK`gX5KP<)ORS zDAjrYV|_(juqcjSzB6ryJ8(aF%blR@H6vck#DwJ9fbo!-f2H$Yy|hpI*xv=IqOO3ami)91fWlO=J>t9 zJ11qXUAx8{GxPq5wUWxo1AiiI&&?!pz&C~1-4NdoJU9W^%*nZ5YWt*l_3DJ3+SfiCr@TZ?mQ@SMWCa>oQfG0KB={@M7jFHD*}U(js77a>Zm!%9*R+CA#9N? zvbP^9Jk&pvZQaJo#}|SXh^_Yv))<=(-q^;>Cm)30eE5X9_wQ>DFFA=bPYuJu!!I?T zQn?qw%ITeX{MNbgQHxF1Rr_J3M4p$;#hw)Mtw}${&D`EcHy@I6j%Zk65(uP^=J6gltz^5sjKCZqVbZ$TAMMrn;b_Tl!(Zob3T+b$ghsiZQm zvVtz!W&Z=dgm>>su=R+ZUJ#fuk>fx~W+N?uls?sZ>T;A?m2ValH$ zOn5A98;wyMdhGvC>z)JI@JPADC(umb7D-Pe6KAaM{l0+^)co|H(wEEt5AUjrCAeGf zZ;$9F8e}tj)3TJ~HsP*%mXsU@j3@fFld_i0Ce%OLe{%B19bey)i~BLtq&og+uM%oj zQpEzD--I&%zCuvtrBOQRBRAzu-u{Hi!N~qCuOC>UFhXysug4k*N$wnjJuMFkW7b!VKlI_e3FS%kdw5<g~HE80qQ03|+%al0BSS)_wi*Awmj)vWzZ#J|Mg*9$M7{L$$#s`SH;qr$>{G_P08>Qa$b&;~J*f6YUnre< z*R`1C>;2QR9I*w>9co7()TVM6{yO}Aa>WKfTb17fJ`Evukm^ypEqg-lvJSAfD)^UY z&~}jY3ZcQIAFJbYaEV2RF6JHNE`4AX{X2W0wePXG$@9$R5;z;B)3B1fB4QJIv)hzk zNuT4uR>_u@osN#>+9#fNI#o-uIjY&(ilLdIJZ@am1jOU#U0@j4bHcoq{FAC~oHn9d z_oQt)DLAmP%K|w|u}VO=)rvb+hH-=+=JucZx{As2*(U0F z>l!*zlQm@+X;JfQS4sAe(iy*>|kkVxYIJfm>fjwi+CRR-Ne_{ZQS1dF;Z%IMbn#-P zUl>yqi!@q3GU7qcIMR|?3ZcP7^_~=>Mnmxmd}h+`%%1IZYe*N3p_DO*&JMj{>a5({ zf#2dzn_+xJNZyNEf!LQ}&Yk)7U9{O=_AO=;3%pc@O9Ko%`@{@df8{O}RCv_sUXu7T z=QI9glF!HgS7Kr!BT5SZ0w`W#&6NDTy7|v7veOXFL_SVyI!a~UK<$l6GTT~PDvDVe za1z0M{x%ztGJxQMh~??0sG5$Vy8=Sx9~h{NViiKx38dTZj+JF8xA{Li(Ir=r2Vn@8 z6n-dJU>zG%zvy1HHGd0x*NEs)ZFbml3-`G~$7h+fFIs!m`6_zV-$>5QaUV`;hR%*( zOsxI`5rfC&V=>HIF+zL z7TJ`u)89_l7q~j#xpN1tAb^b%^3N-GDY1TqNeYIJTCW{XW&`Ci5gthR8F`Q+ZSU8L{ZJpFX7rYxIJ=6cGqmP#9lGtA)YOoTPcU zPV)J6NjeOWc&yv7gW}Q6f_cl1f-Sh&s3$MC(FjOQ{PZY|JmhRM*?W|fnP^g#5&HT{ z<|>~ZX;Wa^HWm3GIEz&oi#l9%OGuA33qKaCzVH3n-QnN*XL_jv zK+-{@?{>k+28wxzXlt(-`&Jr_wB>O!aEg3F&ttxH^1Ge#TaelDAoB@1|x70 zqzQb0obCWP2dFe2k>6kwpP2ZGw0j`9(*fPgnQk~v=rkgBN3AyiT}Dy*3Q;j0Wds7n zT@oaTXV*{5NBAg{7Me+GKN?Trbji;TdSmRvO_($oOA}a}hv#!5%_Jy8 ztmkk4yad?`6-oDd&%LKW2z|gJC@?9jsdL!WRDdwGI2%!R66LJnzzzuKA8lA9^hp>B z+0SWdY9_Ae?hXgy%q*&ALqGPbp=-8fQ3ia7 zG#pcA%e?q9`f;Z@6R>rwq*fb1u~m65P&f0zki?gH%iygK>=MtX z^|p2}K9UiSq_4g|!=`fOmHa*ee!*^R17dqIaFReo(oJa4eGYg1O}(#!1&{5?W}ew% z{pD>6x`goXa11N>IE3|qq|u#tCdmX@#N?~SV?kSFb=udSd-VIS-NMXRwxP#SUgu#k z6*aY%OpPIdw#Sb5}g4H-`b!(lt|DJ~Sk7VYVi! zR4*^R|4~aU*SZF?*UicQVP0E62E{1AcqY9yORC~p`1CSm{NI~{w`9O|Y+6`$6#Qeq z)BSo?9t5ALl+os%X+K+YfBSJc&2O;?gy@D2BaM)-=l3^`U)sv0g}6G?5=Kpx{c* z{mlidfufc-+PgRRCV4j-W$t08CzUNSgB+$x0pK9ea98#q)8ajnbGqr}mZrZE&e2!+ zSIvoHxcw0f=@88{+q-VAq$}B4{uXaB#u;Uc&ZSolG|6F;{;WOpCq5DG`4+j2v zd7k-C7t=G|wqA<3;EN!_4|oU@Q!1qEclyL_${Q1(PvD(qJ~^?2sEe8norRP0MnC9> z*!F*adL)W7i<0rrO@hf|^hZp+Xl}-yIH3ORN9By8l`m-*Fhe9M1XU@WSCM(D{HNvv zJeIIDMl1a=A_9=sa?z6mS8OV_*gJ^RpDG|Y=~(_$aDHwTmMn6wsO$VKN8NN2U#g?_ z%YF7=tv2u7a|T@%`t1Q>Eq2QxM(Jvs&^LxZE;Cm)*_({*eNc>r6Q%w%+iKy%b|fJX zpuos)`7M-3#q`;R6lA?8DLgRrj3V&|2JZ)}Z;ba1(?E56MPgfN@do+jS+ji7QD97f zHdAOUYY^iS%-<|7?6`@+w*%^OF+kz@`1Q#y!y3;0Jlp>3OcLyP!XK0u?>|}BmQwoS zao1NH!za7X030P-YpTCa9q4sS#x-ziGB;r$ihB!KA^l1nsEY+2Iti*88euL%8V`={ z@+IQ2+?#Rv&nhk8aV2>0;AZfVaeuA9Rz0!`NJ_7bvTD!Yj40%ml(dxxeOO$a)WD2l z?+__#j@WSz(o0sNd1M5Qju5tWjegC`?D9X;;cRI@tNt!VdAW3G4Kr$jRoD~x71SIW z%={*|s$dQ$u%^_PsT2iiDst^!!~9KseSPo;JZbvu*l1P8jp#Q=ve?5chhJQc)m-Nt zD;gvCrDdrMpM@C=L5rIbc}fD8_74bC3&-YuWh&m_Jb4z|4rP-n6l~^ATW_dUD>nSO zC2iZi7Ik3UMkoJ~)a{_uB-gREG9EbHIQQ~V&sj=Z4fYABNP18hdg>>YTV=QXem}sV zahzZpE-rG*6Cv(oV8FW8v*yxzOibId`uh7@T@{J73bd1%3YBN9pQ64(-BA_z5N50sxMzSDE-Yd z_1pWVP=uGEatp+CiT$2hO#ZX@dJJ6YimX9$#m%fe+`@O@z=4sL^z?MqSj}B)Sidzi zd;+#6hP*Cgn51LS8e6yKkE3O4?8LEHO>Jljk<(EM)GXvrJ0$50Y35Bo?NbQRNwAxD z>>IfL_2<+NH=}D7%-Dwl6dJLbRe6Kdil2$iTlS_*{SmvAyMv%&3z&;@Ii{<-0ueRE z{3GXTy?|qcR@(==4dm@XfQPHxlW3B0C)NQdg7j)b84Oj1vak~7-qSz@Vy$rbXvH$V z=NO@!6~#H(lPq=j-o1?#_qeajEDzn1nG;RxlGsep4g7)WL=kKjzowe0!_OuKMfM#c zzudony`!V!?EFv)>K^7f(P)~mQLdV%@Y|<)*C9Q8c=$x~ZDOX^orcSEN8YA6fZ1@t z$sPa9z_SeIz3{2)z z6QR5*+3a}_M{{E1;>NhSOhB(fIUA*-!M!7aoeD3Q!!Y8Tu%+?R9A=eqfqv%7QzK0s0|!y|x)>6)U+3k>_z1taRK|ZU;AFKf z(23tq3BFy8&Otm%w!=M>$CGMtbGzAOCV-px4*Nq@%Boa!SV-BLV>C;=GQ% zYy=aQ4~jq4(AN0DVPU_9O7Dx;8)mWOph-Z#p?5(0wFh|U8OhSo@FS68vSb)W^VFCd z^j221v-8>0J)LN+Nzc>?hH3?A$qQ~WQ*^jYF%c~A;X`(mT#D@@S}t2V`7ectZ?oZF zS}MK#*62ETS*4i{MP1%3^6!&hEQGi(QV|u_)>CJnsaPI(_%gBK4n@ino^C^>;;qt8 zo_WVceGJJFVUynwcF1*J95IRztPpFI+uHhh&kqJNpdBAc_6h=6Z@qlVU z=DkGMNN{DWY8jZ(F?M?f`bPNeG_q$w}VMyGH z{5D+9)z7!+|A@VBoS67s=uX4Km0UKphCDk>2e9=B*! zKcjw!Bh|CDMn{Z3^*ZZ^$a&1+NQYHLydRa zpB0IJ-f!eldjKi|spv>{;-0a$W5mmKckCq@oy(VZXliK@+$}EyY#_vS=tw@K`He`& zUOLMbnW%Lfqb4N#atv3vZ;8#9vM8KmU~306_GH2L&_RYD-?xRfIM^=;v<`33HMjhc zOxHU&9at`Uj=l#?#iaE2HTv7&&0N7OH(4n-i!8d8?(hQ^3YB(6dKn&hmrbNw1dKv_ zZvFM|G^FvoRnb3{UU2e5;GU(Z_=87PPf2i{{8$!lHS#l^NNNd(uIsPbqh{QU=kVp9XPfsT_5DJHGs5OFB* z4SaCucVmFMc+^&}OUAhdC<@VZ4GY_sel;{V>^Rwf zY?TFr#c%|oH=zLB;mNpk)FFts_obu)3{&iP4SD@a`ANk*)SOZQ<*XlJ-xI*n-g)Dn zU^}Yd(~Z@c0E~#+2Nd{3xrdUj#t+pko5mIpPMGltfBZOyZjOeU0VtPF;yKEDYTx3+ zOJ9qxAkJY(1V*#P2wkh?xBd7y2t>M4%T^O#7 zgf#v~41!O+9>^UHyxH*s(MtMT$6k*KFB!p=XI|e)9+Ys%vp7|vf7cIF`hmhg_=vd5 zYhTwtu05@I@+8xj7YBVBx(^GnQ%+K|OPYkXI9ncKG zfGWVwFlI4D4~a=ib>0w8TcGB))LM*f*)hp*@lUrdz4(!Dbgw6pe8kTdho6+stl4}m zXgczkrw>YlJVj549~bni(}*AORq^k`bk}51jCCH^xid~l>sI`=r-wR8qkRJUafpW| zMN(j!!XL0>xFdJtaXUVfmE=NJ_Vs%KO`O5@m$^0Ifi@GaEMT#$9$LW`M)f~X;*Qgi zK_Jjn=e>}*(MyMx0j|47d!#8i-)J8xrkR_36ilyt%)T<#G$HG??4L#l5W|s3xn93_ zRfwH20P>X5+zPR4r}uPT*c*!@l>EEQw!^VUMAl0#;pZA=0AtdQliPrGv2k(fkMX$% zu3xKtR>yS4h&$}>)D51tLhiFzl&s4uBkTG3MMdd~x4xpIWAu7?h93bK3FXXMQSRu> z4QoNm;r#@-EEtMKR7OYtTj`dNe^FxpBcr){KPqBZWy5v@xeJIbqVOY@w-GUGLK|M( zaLl{af1dY8sMm(9J$h$Kk9tzI&WNE78eKYt8fQ;8yeATkTC24#`mdu zIyKL~(YC;au|huBB;3{IY}mj~zscgAZHdj=E~^$r$3K12@11FHFsZZgxQNv#ZKd-# ztk^YoK$NFFuo zkA{|zdJssFT8VuV^&jA*{9`nGrMw59YMMA_^GN!zn5^4-7X`=nm{#*+M1g`e zGmSUidtx0lDmqCWe2%AZIa${`f)^KE&b$6?GIZT4 z(rJ3TwCcIp)sEB++A-x_pWn(Yezc>y;}Lk=fu5dvl2JW%d`VZVzvJ8cT~k-sd;3lF z$}1ipm|C~E*cKtQONK?M-tziB$5({Rb+cBY5;n_-=Bsp3m0kdCW!NOZg9z6LKre!V zA9m*VbyxRZ3ogt^J51b_-_bNIy8ju?8KF*Oe-4pdl3lUI`(1`dl8*x6S(Z26V0=a| zK>k=EZFp<})!srr^(Wh%eh#t+@8_*df7%z2gNCDUagrgqVBVaOTauA}-$z!?wUn11 zu=Stn6Wk`w7-4j9(?=&GPf)0(b1dNf8A84rH+3_0vF5af&3+b^S_kSXvvNUkIi8sl zk-i3FPyK84EJAagHj6xr*y{0ZkA4!F^n#KFd!ZuilQVjHyH?Ti#z9xKvHwvT2=Le2 zcIv1e8xiuIERNq(3X)}$6Wi{MUT8eX$PHkh#N$`MU|zhqhp0tjE(h*W(nd>cZ+JTp z-X6F%dIySBHAO4_83DR_p|&-ZX-!r7CgXMgjwsQ+xMUXqELz`AhFp7fA+8;ueo<0( z`(V=^EySdE?~e6GO4EbhH2jmh&9ujME}Y8tP3N}Y(GI5cD*s{Sg)j^00MF*_U1-dvxc0?0nTI`n zYR+VM#7-dI8teb|SOR-*>IJ4xx38Uq6RoJ3()^a6*1>7~xU}t^U$d_p5lRDL!F~J} zT*yxcK4GP75UcNTQpEnS=-&yQx9s0BPYo`3AYCucM%aLkmXWT zjACw!4nX0)$Bqy&7oaa<75?p{R8>^V5X!l{r&KCQJ{44`ln9Ck?u?;q_!?ANrm2!HH8}j?LyV1!5&{DyRbc0Td#wE)H0+%K^H;6lqjlaw zl;Y5Rcjv0bK78fDZwskug`7a_QuRZcCaZXhV?H~NQnU1FiTLYNg{V`Ib_eVeRy@7w z7KS5Atj3Z?pRJ{D+wcYR$@<#u(dD0t8=GjY?)x#nuwzp~La4g{t^4&0A+x*_1Gd-r z7BhD-;M6y2=8-nf?vkk=>8{*vOI@kulIg*D1cpZo9N-*({j8v z(zga!@kIP!JnLXzq$D20ufRh4Z@7qWNXM+pImNU-|KL~G*&PM50}4M1Eoh@n!-e*} zy;HR*Ag@WQKG(uNabA*YyK28}O%AQi)eiZ^4~x#z9~S9W zbRn;(9XaB!d*ze33J`c(?Ccnzoq*U^Z@0jtGVm$VSaL|h(>y!*onoWD!*swGOG!a! z0}0s^&2{tl>tIWqp`$;a@L}_H*AMI~3vP$xmrSVSFDx&Vjres;RM|T%&eWWBnEj*6 zx_P5FqLPEtKXQ*BRoG2@rv$7sYs|1urLbZQ)suCiSP{w z3i{NfYs93%Zj0W9<9PontvYF~($()bMO!!SN=r)%4G;GLH!zNR$KAznbt-1?7mDz< zwvaAHW5o$Ss_xyLUUlc`;BM`+hKhD)FC`_tOIu(5lndPiHY54I?%9&7iJny#RZgCi zc$ctSL+irGMBq)`jhvi4#0mmUq(6#7PQW@UZf^4Jt}_?p{{1kU#k365c$1mOWDuMW3JCwb`$Doq>5rvMwZXioNlBo|>b0FC)v=#ICNcW3v+v0n;TTPw8K z)2%HcB)!S8g6Z5`(~9Vg0g8C;ob#c@#zmyP$54k0tS|0Jm?ppxi+oIcv{mDMoShMa z&g*taYx)%yB0VL9To8M~wu`w?n75Igy%S_T!bt_TkguOou>~wXaESNi&vpl%%e}Fk zh!&Ln^ACgAs3meYz9}r&8G$mNVOSXrG<3HZI17ylun^${ln(Vssul*o>1@IAYrSXT z4}UBgmGkWIBC|-&rCIeV#Y`hvQIj@=hR66V&CpxfCPwMTlkBQsIFNBU?wC5JuH&;k zgtf>^K0?9LKLbQQ4#J+pHsz=D`fH!*x^P<#tWlEUM8QbDHTp|$dD>bqK!A9Bud)_V z-`ap8>U2>?^*AjLE(^UeKunI6#i9|keu#lYAl0J^zO z)MXWbNJNc`Ci579di(lJpvVBx{~0)psPF4mKfFMv=;hCw z95=u;%NA2f$D8TIA1fd*kPY-Abk&abfrD*%QiPdN?@@{onCHsUn7acOu^7tr?2gRK zMa8MG`y$+RxF1r%zl!KWpz||>N#A2hD=t(46KE?s)9d7wF=`m9IForXu6IDB-jAusM0O_d% zM(!o8tG^2~6?1LP2FN5CX$bkIx3>23++-hs^PO)_M+Jk7UYYlio{HIS>T~#uG6+m3 zK+giV>{z=huPLx3Jc1rV_sy~GdG9b`N(jb~wIBneo_c~wIU*I|l@eM4!VVF=1aj{4 zEjGsCs8_#gKZfqWvL$03!Nx`6fbo4#GME6J784xaqK&%eyHbsnk263bgH;OLvoUew z{r{vV`T^=GqgZByU@|-FQ1~%)TJnG9X&45D@E-aK0XnSRsIXcN+t|Ff`^yMy1{;JZ zI)I#lQCx&s?-?B(Uku}$^qfmAC7VB7)e&N+Ter@3ev>*q0|QPaAgsS%j{Pr~A_KaP zV>^%$xDb1QeoZLO(8Jf9V2BzemOA+nP*|D^{+c~!b#;4?qcNt#({>0rX@+q z5{wotEvodNmTfY2gLgn>{&6pTsSPunVT_)kLprHHuf?uP=Dj4iVZURtyy(QYj+vRT zDp#p^2)_TNBAi9jQ_sb0`q4 ztnUF(6XyanND8ALXAH?XPAMP%aS`%C?}O8(`2TY_=Z3IIiE|C1Y@9iNJ_yfw3t>G2 zzP5$xOKxE}J?HOt16#~PUt#%CcpzqU6KMg7MSzYhE(BxhZGwT90|^JjP}1K+edr{x991QcH#Q?NKv6!^%!cd9OeSPPWCU ze^%`!kiUs^6i+b}%AdhObHooEc{4rK%z|moLBh5Gz|n&wfFoNzJr*;2f2EH1e`~Q(BDn|xaZU`VK9RvAcBUB^ zlXNUB71%hyOuXS4B>+VYhRVddkSNSRw}Bi8pOt6=f@McZ5R7q;0|lB;SqOLtj9s5a zXgLw0N2ibNT#H+A4bULl^S;pRx;FY$v;GG||C{yf^mvU62Fr@&`PuI+TBH0==B#aI z8_2n3At~yGsN$yQGun9aOTh)Z^&gr z1_GRc$_S&PL#yYdQBodeuoPgY2a3Sz=}FkQl?Rv%RAlNaf>fM7@eJlp+2M1 ztDMK{>pv643PJ#A;xba#&Jt*>aWK6NfshEZN!mbJcRdR+C?PSC@SVUaU38LQCZ9?` z;I8@s9{g2efJJDgF>3d2I1jqO5@xLug$e{wy~d?2;dp;V`PZB2C{)-ZWxWL#?0YPW zcjbA7QF?kF8MZsObG z#}t7em6Qm6O&ALw_}PDqePeHz=R~%IMNt=4jsxf&6cqThP&?jdNtLl5<#Ah_P)8Y& z_2uanvOsI4I)V_Uhv+$LW4BW=qo$Sy(1iHygmko<{t}{T*oqkn11Mt)%3~<~K*+?f zfC>p6UKamN2=2dG%i|G9%}qr`#kQU+z@zB*$60EMva)Rq!#@^xa^y>K8@L-W_y|w{ zmx*rf)9w7g79;8~dqilB_xMl)V?51Yi!3hEv5{m<)N+ecQ7e#@K4eL~i7dWafY!M$ zdjqkkⅈ6DHDP>{8{t#{hD)Qk%kXT45-U3^JH60a>d!ki z?5?iX(AWNzXrO1i?pWY)4NdK!$%gg*mk+dLjCd?JGnPR+Ue?HIRGM>_r6f+s59NSxUxaWUwPwe)k4W!URggmeHNK zMbfA~5KsItdLZVv7#Tv^@c#YRlL|YWHrDGVaq|DXqM$AKQ&2NYpGuWJ2(>)YR0!1yYz|mjF|^ zfS)6;06HO<3KfeS1xEQv2H9jFL}p!ZD!D4vW(vhLfyp3N14)Reb-Ib~DWR}~@xz1B zJgjXoVRcxXs1Um~a3*hT*R{8d8Vat0_-nj1Fhla?(l-%swGQ3k9iB*mZ+<) zNxWO>SfO9wYY(Z_D-_tj8@RqWNA%TlpMV<~79DIRY;2+}7wOkh(@I$sA=_ZWMnk-J zVeR9I*!n@+xwOfnl!Q852<^&3#|ohrbt$aOv5=5MSY@QkBhDXal-iPu2Jm2sER<*s zb0*5U#jIMyiAxgpZdrBs9%Q%o&NJ3z@orLZ%m3yTB)rjvu!`DUrcNT^b)o3WT=0O# zQ(*EEl#t};52+sP*~uL@9-uhf->E1cqql~csOrG0fc($En>m?=N2L7)XQ=(6&?Djb zO+HsnDccRb0NNMg?QDN|tB_3vl@=la#bY^XJ0;e&y}x-%XN;CvxW6Z=OvcN;a_gqf zuN-u2s5T&y-^m;v^RKiUljs8NXnp9ZF$ zpiG)#bWBpO4VKV%{d&}Gw1}`J#6cFR&Y<<=jtvABAhD8E$rmraZtJnk)Gk^*1`VwU zF()G-WYwaHq_Iiv|F!m*WCgqVS?<_iA#)-0Kty^4T8?3R(}}hIthS3kOomsSxz*a` zDeZqndDU;wv`HG|P^s8?n2nY0+?%s>UL}($Sgw>x@7sF|k#>oJ zkkTP%Y|3SPBltxi6p@$&142YcEwdhn3~oh3HIjDWb&CDyue)TqkhbP`UDS2X*X!3( zBJM$joMR_z+kGc<_H~|@;3%kq5nX6%K4BjA%yripvdsmeA2mRujtb`$=IU@lt1W0L zQ5B$|w7vV$rc9!_tADJp-1*vp_pT3k&xUh;El|FK)k?XsFLu_fZTEn@%Z7wUdrLU^ z3cc+Ooz{1_PbjcQ3gcLB)-O>9efT)`eR4N+ThI5!jafl>VU17!7e%(0iq`)$T$T!F zCcYd$-5n!kR5E6v6CJ(SePO0_;kMFD@sPtp?9Np4s(zt8B6m8=|7<^b7McUlXJOG! zhYDK_#Za&EVBg+OqI~%LOo|1ktB`t+=zYOv$rHnPW5+p#0Er9%|3P)gI-5YGswuoaW^jWf8vN; z0->W3wK=JC37L}Y6cTLe(Izg!RmlNnbU3q3^inlDP;`hyx-m8;N?k(JAaW_9U;wAn z!qp808piU>2M-=3?$$(PWD#}#X0r{f;^7T|BWWr`t$#_!mF*>+>f&9YfcW^F$6JT* zLBUUhB=>iy8Ac;IDk^tyJaMY8_sZ2r2E!odNn{|hpqum*iy^qga{k2l}=kdP= zgCjY;`z`efIgCEL(_-JMnWog(Gd}|6aGl`wLeA{A6eYx)+txolKDK5X&lo+! zp~7Q7^4*00jkO2g?>^icCE8h(Y&Gj#!}p^=i~)s`pBe4UtDJb&gB8ps6;XOZ^6t!1 z!5_J$+T2)@bLBU;nRbhO_Z{~te1%GEadD~+qKlkAZ=S94Dnp&@#@e^k$xSgrF+1F1 zVQvl%kT=pl(JA3l08^+!Q8)5N5fwVIsU{~OptELV`~-pnf!u~}JjV7%6Iq63XW2%$ zA=4wRoM`e4?Gwq z@@~+(qG3lVzhms5(4UfjD*sFt3x2C-jQ!YrH2arlMWp-2?6$b#t|fP~<}3xa&23D! zhn!xS^&h!4!h5eX)EMm=|9%+`R8IYM951z+kHAyLHRnpT_6euOgK{h7azT-z>jFF| zkOjK4oV7%i4s2wpJ5rXO;^uGR*KsGU!4&mfiC{~y<$f|;vTC)zQ^N3u0?k(9`P}s+ z%)w?3;YEV;&j+u@(vO*r)(7;^X5L(PuzM%z|jF#A?%3LX93r zkzJMWJaH@lo&$lSGXJUD*=4p(Q7Phvbb4OWeTHN}Dk@S7TIY1pfN%QWkQX3D=H-uv zYwGGM0LiDCH@)M??QtHxCBHnc6_Oq|MX8~pZMlZVFz)dA5JMpm8J6^;DktqWifBLV zIjnk6r9NnKQuATNqiENPQv*Dk?+(7(opBB0jg&1xO4$0BmQvIrtlLmL1|U_bsM%E% z4JDTr4c)gnI_7q@wLlOcpoWfSkk{e{Q7#KoB?hEvK77RfD&sY4^pC|cAClLJppaxg zlTzQ@k6$2m6~bQe#p!+ov>?`f@OQlnbc`@$3zz5jpd3W+ivqPvI{u^$s z-Y~RP(B=|JJNR`#%Jm$x)sNO_{{GEVMmSuh8?p*A$zqAg(@UU@7t8UwD|h!CTHTx-<^pUQrD}qj20(bxtrDn3lDtF-GdjH zn6mHauIq3s?Q+c{$2Cd24YQhVaPB_bM?K)aOi6ot#V3aLAAkRV((tl{Alc?z$2C6v zCWa$G7ujG8SwX;W3}LfVrXd_Bw~ycfi1(|8gJduOTZ9~toiKg`kanloo0c$>NjN0_ zqkD}Ccz}?S0t+5L7`gEg2LJ*@5m1P5mI=*u%zY&g2*@hC_~Q)k?cpC=-PZ{dRFahJ zWQ3FVl!HJY@Cs5*XTSFlC;@DTdSIMW3gX&0B+lZ04!1T%hk5Z+2-xhG;?B#qQn}Pe zZTUr?3i^^a7i7;%nTq>%xX7N^A?Rk1$!MYL5;*QJDf(SoP*iAbbB8pSR7f7e1pP%m zdH+j+^5u8U)LFTBgRSdRb8Ls#PIF7I#Y>jVM7akw+&4uPuKaWA#xJKC?}1;L+~fuI z6!xhrBKx#t_y~xyP8<>`tF`|XSQ&pV?k}yJp7iX@9`6COA9|*Nkq1gh% z28V~ZoZz=eK1N{Z?clxoCOQneAIL3GQ@rvy?+A|Vq>`x{{J-rs_u=Re#*c<$(_D95 z)YAdz=O)PGG!?vJConAYM;_jvEcX|qM||^v{>rJSwQDI~y)qeSjO8v7XnC0y%{ZhZ zQuh0=h3MjMS(|Hxl)z@UqCMa+U>*Q&t5ExIqvU#4}~F-)&m(tWpLgb zH~%%hBxoo=m|wo z)|3y|zqvp=J$!S!lVLkw>beO2ZQ|_ZPtAd)1TOyKdC0yIb0#&!uk|9ij9p!9xPio= zk@8k>P|5a>3o*J;d)N|&?{cW1{)2R8Co}sdI8jnl*QZr`2^sE{+4Xm`uoC?_5w(Pb-TVHXJ^Sj- z8?GHaWkGiz2?(FA%5a}S?NN?RZ-Xgz5SIi!@-0Z{t6xcKfoplTktqXllfcTX$$J|9 z)6)0w*3N3BQc<(-?qX9fPK#Y@nId(@CJXc4(z`quNF0Xi4bdB7T>a3Khpa{c^Co@q z&ab&QuOyyJT|*8H!NjAlt@ebX%Cp-r$MdTkkiPYY!{_VxILEV5#R64!RSb$Yv&fB8 zpbNKEQCI(u<^^}c-NIZp#NrGU*EkA*)dNLGFJF;-hKs(3`5`_AXsRexhq3(N0yJd4 z8Ha~aV`TqKn=$;q-i9w@;H(COhQvJfm&;Tw;(oaO2FH(pVc{-dI{GUYZLm=~=i$|@ zBrvIi_oVPZgGD-NVw{UIf=lD)O7}zB`qNVWoMX`j1p@mt5mEs|fXGgaexTX>KhHN4 zR9^}T6y_KZ6oY?_BVz|wEXL}3w!a-*m{@5_uM=gAFyS_SY``S=jDb_{VP=B(wt(Wz z%^?D>6N;|d+#PpRTzjJ7tlNE~6XJZ{T=y$ygN}2t&$z!YaoeQk72|uYH zZho}+iOCx`_Y)s$zQrXt<+aV+?pSb==UsZ78}}~3VkCJg`DQyU%C8{#<8OF8&WkMD zGZ%zD7ztca{My)@I<9-(Wf^quAJL4RWulqzR~L+}V?yvnJxg z?alJzo-Cci>2>88j+GkU@}At+^V#<3B!0c}r;3=X5cdEwe!_ZeiA*4iP@F4RZt@5bbCENl z1e`+Z02)IrrK&kq?l`%m35@h1j|;Lw^?+}u-((NPdw+L$Rs#l+ZH-hS{n zDnP36^0idVF?DHs@wc5mCZCQLJD^EPOEKG&;^uSOeg4}I%*e-gx=-$e-Or(Ez27fR zK8q}qDwsIA*h#C&|$Fgc7oe+9OI&`li3J{<=QbeAw8pvS+BIkVK)G>BdGH83Yk*Nr1 znKSb*)Vvs$Osp;VmzmP#%WX}3zdEIvA(DJscZ}Ri!5~V%`U;QIUt%{XX`6rZe z0zUgvKdJjb)Er~M-=ynZ&4tSnQRCl&TJXULSquQm>g$hWd>aR7992xS8V#bE)nq!! zWa+u%?fSA&95-@#@LfNKefg|oF8#8pq#uTMoUW7=1BK=*JD+C3&8O#!?ARyHf^ll{ zBkA5yB(+R8Cf_tB~+uOikAS2CPf98`9wcYDq!Snl%9>?$GrZwsf4t z*9FpkV<1`0>lvkkEvsgvvmN5rEd2Y+ra@_uNPqvuIf$*);$V9`V2ItkA73aQx;oN5! zgsn!6#U$3g{4d+=XWAu#{tUNba zXl!hMP4-eTQ;@U#PMkVblD2mko+TuolEu;)uqU|)i2F@)q}IHFKQIRqXn-6btFf&S z&Nc=VOgwMLY2C+Z|98PO`lT{S$KzS!L>@03a@T|MyrN%W#c>g*m$RTxz&4rsMS=u5 z_lkHQMGKX{W^*FgKA+0{-xYUXY%h%{a?}aF=EWZ!=NEHnqfWadK_YWv(5 zdsFnV4tTz)I2?RQSUiNdF_0m?Etm_P<=odlVMy&y)*JMm^;b?xBrZ}EjKyfa2X{-*dCK}aunj*%c}aP7Bxr0@L4eAMNjJc zW>WKQ<20en)LzHBisYv=)_Y|5>-Bk&(k!cwkJWxXAJBZI4+O>V?>pVCQFS1F5RGLH zg41|zohviuREJ3Aj~hg*s&S|GnN$xrF69VJggZ~x|IX(&Ylv1ff9l29u*?t96DFGk z7)u}RJKdTEuT3xhtGgkQ77%!9@jCW8N2#d(Xap`l*TdVqK3bK;l?mK2Zfq2|$p2b} zVz5lelfiuVZ}!*6fzhi8Z!~}El*@>PQK@dh~dIkB4hYV+AmaNj> zXuhS`b-K&7!_Q=s=*jfA&&eXVDzZDG#!2IFq3)Xru{bJbMEW91S{{#?{Qg5Z1JQxw zS=k?%uF~RxH+)KoM6ao1I(uldtT$c1DbvYM7BxP9xI@2e z$6e%7UfJ8jgkB!B&wZ|D7LGh2Deg-{4zpc;f3xd1g<2o@>$HX38IlHI{bmNCkn~Vu zIugdCeMZTJR@;u)(|;4kU6!k(>K?Gp)IYgYvRV#3O=#e4Nwq6=c*aRddsZZ(G+H%K z^x(m}m`rcP>gFj}8oAy>A}*E|*-QulVan#gIS8SGjLM76?n|5n)3N)yU|3a#9pC{h zv!`v$O-iI-)-->iZLB=OXvv%ihus&BgpRV`r-;+uU7Te_!?k+D4I2y&zPBr#eP)j= zd@jdk>yDi}OI9m+b)Qfh@$2h zcQrCt!}OGgw!+p}Pv^#oZM1$o8=gpw#0HBAmxvCV%VQfv=bqAEOtUNfC7MFTED%v+ zQtU-daLDf|Mj`h3JzGecpFlnoU{ofJmO++f@u8Vt-h9jBeuNq{t@<^TkP>omme<6L z`^wB}I>QI2pt#-PYXAJ^mxDQ-UJT>NVVmb%+%D#w8hG6hkQ4nxayqu(5z319R@`r{ zm1J452bNuI3=n11*vsK5Evk{EGv_0xQZOcnTdMgM!?>5)yI<)i?Wib9+uJ)e_lDo| z+EW+-9QDBY(HY_!*!K(E?RTuo4U@znZ*2AFx6hOnjG396)1wxt@_wALaEXg0H~v{x z$J7aSRdSl5MoKfu*hmtlFnE~70@M9p=vm2R_O<+pRM@m6b zgMrG$t}s`1dc4MHc|CUgH%|2}#9d{!i^bhP#qDpx>~x&MJ@6^khsaIRvAU-sc8wL6 zaT)~y^0}_##k6#lbMJsk*1yz zfe|E2{#@<#a$+GCg0M^KqxbNNGs9RR2gXDiA`!M^U;_)1;uYj;V;na+FF;oWx{p6D zm<0ycq*)|kmTaw#rBVA;+ev7zDu``4{tgv|;nJXPPtbucMG$iN-Mzcc)%K?5t20;U zTSeh|8OpffwNz_d7mB{*p*5S8#s4stUOjkQQ>5)eBM-|52*kwVNz@e|_2NIWhH1U_ zK%=?RXxzJ|-B7f6@f$OzixuNaIQE|Nz!3#q?@O6Wcf~#T0;3B0WWGQS;S%|UVQ#Q50-*?NsAy*bLSHMz}*7-qX@UN56jC_{W&eCe< z$mV~Z%(e2UG(1REP8vi)UnO-HAq$h(LKLjTm%wd0hHo{!_U$pci}Pi><70wLFqiNn zOOR;Z;G2f~e)ph8t{ZA2JO%qs(lLF)IOrNC_$q>jImMxD4=JU4a29XV&U#iR89(n{YozGy`@Aa$cZ7YHIDF&a zYhko-h)bV6NIYj^zA62gxo!%V`mV{m##OWM;*%FQtm=nN82l)4%HA$!xpBo4Cr@@g z)H+qHp}xh3C~2cy{~=JV7VZE!5zM->3zvimysQ@Exw?;j^~?LYFZa8b6Jtp-mn6kA zj4bZK{fGkGc~Q~^pvKZ#snp9v14hER;sR2P0?MfZYUBw$u*<0SnbM4-(;Q{Fy=OsN8#bkEJK&#+%`1+v^k%K zEcw+69a-mG##J^rqTH*hIywwQm-MZqRU;)qfBpZEcF=?PG#v=yAq}5Je}5IC9o`BdeO-$IqXCQkS1(<8c~zMfiKR zn9|`5cy3f~t+GI$JFgxRZfGFw0wz-h=o4VTSeI|dnBXm zjBZ;-_6j93QYp8Yoh^jQ%FLF1CzR|hO0su$Rzi{;!gF3ezwh^Xp5yudaeR)?(NWy) z`+8s3d7bC$e7#<0+$@tI86lV{Kt8nX4TG1nB!fn-v4(yiB!;+9J8Vnie_iP3If_NW z_)U-6Wto>gfj_Bi)U4+qTY{<+rcEd+3J3BrKDbF6(J*Ms6R1O%ka2opchQ}>?(Zi7 zlqUvxUM0-K8~au|>_7TNnf2doEz-5PKK8RgPNDr{&;7Do{$FK?T#53;4t* z&Z$H;hOWzIc_*~FzuK}8XzOWfX-y}jIK_BXir`wcrB^vV@lNAYp*PoZ&!gMW@*SWb zhEeevE4pubSFqI_m5g}b!PoB7+FaWz6xPaVVGpRQ^~nrrDaZXQblG6L2BCQb;B);L z1s-&*BRf2p)qYg0f=s*V@L=+g;Rw$H$;(}&bAUsRLKlEG=O8a<-va#~FN|uzL)2|3 z=;|BcMDrkf=%dsQMLs+RrYi zf7|W3g0}hfyqo#&c=P9K<19+UfF2}Q{^)8=z4xRtZIErzxZQ0r);Z7DegDwKbjxN+ z&`F}yPG!mNW!s#QQM{Vaiq^I!+z5!Vc(9OmZ~}#$jaCs-h+!)IapZ3n*o6(>-NA=E z5c3Oq?_%Hm2;V^-f$Cues9AvIfUF!*pb&$U_r2YC)%_80pALt#423qZUlKzTQov;x z9}i_Gv_UUtjc(O?xbzB{ql^m{0V*z`q5$$tR-Mm0apr*@IF(=rUHKFGc?FE}6fLh* zPf^o33=yA@(FmpwTPYo8g^OpiYh4P4{z2l{!tBYrAfZ|teN~<->H_bJd3Q}dXX!5W zDiYGkZE7Mpg93S`mPVviu=m{?MX<*vd^*}kwElu*)paD;-*opMafT`|u)Dym=Wnlnh)=S5F@S+N~ z;-g?Q%7>0XRJ2bJlMeFeM~uV|Q(Vri4DSIrz}seZd;TS{=rXfSY};JicYhyW_U*Uh z?BA|vo(2b0k&mIT+ps?X6(kCuxu5zyV3AZSvp~$}!8OBPenqwv=hq}NGWPo=R#dw*m}#yh^JMZ0lm3L_gg!^pj+n|bKsBC+s-;_4!kx4nnA5dI+R877#w&2E&cemBPGpB9+=>=(jyiCFa+Gr(; zCo6kU4Uq_cj+~3`SnX5q{Q_G`iq1TXCCab6Y`JjvjSJiUSr?&S_Xc-yALE znArfT;$nFWTwP>&8jtVSiNv>UdTH~zZ)o&2Pk63`QTqOE+T|4qfZfA}6Z?AYLL>3M zd`8a@FmU}z`5TC3|H^smDw`+chD18BF0ZNTR+=0b_45bfSbCut71L6mcp|9vIPs<( z?QyS3m_n!$#s8FP_eV#}TfYF3XVPq8hZlz3}WeGYPf#mCx7 zo9{G#-)oz5)TS&0&@jbq0%eB;c>GmA*&x`+LqI2ifc|pUT4&EYKa-FOcbv&4@~R(E z#Rh9UHWAlzw`*OuzJG2XfP4UeNoM|0W>VRFiE5hzfE>m)em*7-)$cHPMNxzS4t(?4d z20EzRv1)BzRlt-#^9Uh@NjZJ(cN@%AW~c8y4Xv7LnIScRC&1dE5%oPVx!&rX#sN#v zdDtIawh-y*-LVDYF3OA!L@U<2G3P9Ga!W5<*rN;TaZK*N`gE!4%bepQ#@spKUMVPo zHSvQXDvN05CVcX1id10_@9Evandk|pD7u*pr3~B&xES9Cu|bZiVS7L49A${FRMk!~ z!qpOOv9X2-7(q-_bamuvXHbkBqYsp7R#X4Zx1IP?=puW|dp#-zqicEVFzHM`A>aW$ z8~|0BG~#JaX^{gDReJYkYUS_uSoUZSX0;C#(_++i3&`2*o6=RcpE*RviJT(*m4VDd zXVj~nq5AvCC!df)3=iGbkJAM28*W#cOG$aYztvtd1hRf?R*Bu27uOJ0thpV)vfCZ} z4u*<|bSMi4^Dzdoh6-mxFH>`*^5-6?WQP{s346&_P9v*9)G9 zqVZqy9(eE&V=s^XHvKU4&rOc%JJiTR;S!5^SvG)&kzHvA8ClTPDm^NfX?JFgzF)g9 zVvbNKV-@Ci3A8z5E+#0q)RYT_%r_6qoahd_=ZWGUTG5pHss)_LVgXDAgRh)?WBv;y z{;!WSW$k@}s<2nz=bL%3=)NY(2=Svi6_yi;d&UNB@{Dp0%tyvw`L$(C76 zZbU@{Q66-+Hn@P=-GL&z-J}Q1gcO9Irv!1g>sFl3APEs$;=0XfknJcucI=n7g{V7C z{rzWha&i)LYC${E0vYn4{!s#o4@Pu&)q3roDOLJkBBjYR(CO{*(&0S^FxttdAFoy! z{H;qPVgL;In+VcNdcDKIyeL}1A(z{UYg7yW2NCP>V!fb;a1e-yD}xkJ@;fj>)P!?g zM#O;k9hwgBv1~cc?hM|e5p3>M$c6Ck7uY!HiCq<95w4-VZ`Z_Hrr$kjN!Y)yHc7#_$ zWpZcqit*4V!-pxc^p+^KoPeHz;5R{kShq+)GU}l7JW>7c)vr+zT|D=~G5GPL!5J=>i( z=72|WyAi!x>gA}ec|hA%m3=pq9|My8;uz_Oft~?WT+p(h28jhTS|q?F7Cv-`So*T# ztfPKHy+UGJr4MScR2bq5UzY&q?-ZLML-pSOAmGRz1NAh}k3;{K6x4b$GP^wy(o2J3-~LM!3siI0BHPZe=v!|o#(38P^&3qL6y9X~iMk4)lY zRji`onx5?ZceX5L%p1%Tzoqy)8=7!?O+GyMuf>>%fhvNNjw1T3s+F?RW%K{}YXS7v zliTLB$KjDEuIM)YBmeKWPQq^`N9_d#*W7eYtO@t4`lC}A@&hI>=Ezm)0lv5x9g+z1 zQS2~BKVotx|9fQselRxDXUEG?GX*}v|M?e@qFJ`u58^HyWC{R-!Lzk0eWRTt{<+K+ zHL}6DDO(2j3C<^DOJj@7)@9bt^(nYzO@k}hNC@D01e!x0m$-HV(nQC~l@SW2UaRY-ptZNJm6bP(h{)& z0Fq1r4x9I>FuLQ6yHFrNF1S?nw!59MmWUCRgx|`>&pd*B#YsqTP7F-+uT~n{;l7E7 z6T*8aq&!jAzqFZ?62l+-;of@A@)vWb_;#M+sD>GtFr z6&!yjq+kM8J5I$ciMtBv#0*rjCWGrYp%!J$Hk&Xe{}P(-SQ&cp`R6UsA$zUQ`+xty z{XMhY>s$UTqzb^mG~=PUcma#&F+mDWm(8Ts`) zGRyfq?Zf4ufK6`}SfaG9rsB78M1fI)r|_h@lAKickACo{8w}N>FRnae7*r*|Y*@4xMBG5zBSI9_1U42Om1?E4nUT#n3GAot!efyWrT#}R z_JbK_Hy+nm)`jfeO5M-l!RQsTCcl+EK2_Ex!N8S0TqBrJ`&C#djyR)u&X#x^{qhV! zGCe1zP)WLzed|(0a(<{gAj8;#KYjkPOYa-Gl?O+iAVsIs0759RX&=s84LyP7b*voI zjT1F3$0pha_pqav%cz7=ei`_Xg@t6~;Yb^eA0$2yvQn?MwZ(+}-|sIdtH-h(X09D= zaK#ju(I&z!WOOuq7A}cqKf)d=mH3l%@=-@WVp-RE?fRL8r;7dU|F5l$>WQNk7$8-E z+=>&;lD)TO_+MbW^ZMl~uvy`4qu*PqCizf44kjMDWbei_c$ixv9?) zGs5HY^+^a#R1vz@L{9Uz$3ci7sH^Pl^gy3Bs8y#G8W1G@PsTmrw);yrd=8u^VT_YVEF?`secXm|kh9;+?<*ww zs+Z$YD^nn|3;oz0OLp1jPj@dHV!VxB>-fR5_lLm5XBcEo3xT}V0}d$Q^fg8L+@)_b z;ROb)Nc>rMq?!bwA9kWy`?*W%JO_q)Emi(?w@A_3B#(P9w|0to6Zn^sb962?2!>)= zJ}fMjK5ar#vB-PbnzpCSTI()p@B-B5J)h4}2)8&evNA8ajUgs5YSTX`c^dF9yt*(c zhN;5!RD|@1%7H0H{FxJD?n3I*&&0RIx;50Tr~cryy-&)5Lub;DJ`fK%Q~->+c{fqB zhs4nYBSe;odkgN-g-c2~SS5h7utu1#bB;u;Z*J!AW~B}+ZRQ#CFMaZMFu@Q|RTk!J zKhH2J9f+pe_E zL(s^B77zS=$=h`e?sG3*&WL$9yqO7?8WI?W#4U{ zG-lo(B23CW^tj{XysYV_JT&3#w(+dKQp%%kbdrYYovVp+cU7ynn{6jCyq&R#Q{>Uu zM85ekKxW4WUMRJFDSZUlLg7+Q8FOUr%^WQFd|uI9QcN8q&v@Fz+Uj8;{R~bGZq(ZF z96G8SOPTokJ#~(R8<+{@9@%*Cpl+CD?w* zU&o)LfYbz`deGu|2R@#FUBiE@GO31#+B{2}hggg8f3gK{lMZuZ7Y^$|HoW@Ppzc%? zY=FL&a}=3j04!467IQuqbpB+{vtS1DsCN57RNWUqTozYlh3!qoW-cpumciep(x z9?WA>Jm9B9`bGxiENepb8|3Un{BON|*>8V1~kX~=aN#F6` zn6C9Gh7G%k3uedX94oQ)Kh}3<)#d9G{zw&Nt6Eh15yrT*Cd$~rec{jA0B{@Cwd=jv z@`t`ghJ%NFUyBdYV3zq|NwMi++RdE9LrDD|HkfNT@5<2Akco9WU;f^m>SgfD_q&AG z)SJVG1LDHIhuEcOEzF;lJAFvc2^jg&2T}_e3!wdlb-?YdXb55%+`xh3z{8Ca4%IiUT+a5zvc)X-f~rEp z^~4;yO^*uO(Qv$y<*%$Q8k3WTExJ{Af8Lk<*+$)~C}F)9R1nX9Rz&h!_T60xZ%Ii})ha7ev*reYQ+%`ZQ9%<-W9dvqb;7^$b2lDfZSg0b=S- z<^8r%29khLi*&G)mW)WqaV!n}J|RCc|Mg0#XBvBA|F+8|t#df8fVcYG-qHzr5t7s| zzlwCsR^LASt$+A?zB<^Faeu44{u6Y~Uez4O(V8inRV>PGs7enXHn{w4&OJB*3hyJD zp2ml_w|L9Eau>IsI|^uukYtO!^cq}wOdE7zKrN*T6NMk_de)da$Xv+d3D6;+!X6Yb z^INo}kVs%LvKZdpFT-T#U^!!LO379wAlrvzR{tME*rbJ}U9K*!VZA2KfUB51e1=s= zDN@H}q^OttLRuGBM|`)%oQ!os$hLYqdGJ}dCUvqH;K+J`gjVhZ+e7Sjt(M|ABte)V zt3HY=A^ViZ!A`OsC+hZU51Y|r(atese=^3Kl5=`fZPD2 z2|q-O{myxBQEG#KlH(OQs(Tvl(}Q7m-$R-;wzk+TNLEq!;+^bBMHZ7!qu@4HZS|6LJ#h^uKGK45P;;h}#_$&rb8%`?8pN{%hUM(l0;! zc@q$T1$)WU#D$WyhyaDmSG|QgGF&Hmv`Z3WXENCwB6wzs#j)K5bz)>Khp^7U<-he#e?o$+xuk6mv<0 z+h1%rV>kU?@|dF&xn4_37ENx|zH<5W*RI1b)3LjUghNxWPM>`Idxz#vhrQ2f)ct+< z+AwP1`T6{0%c^1{@js1?j~Q@&0b{&nZJw;+_{Da-o*ZV0fD<;Uc3yLA8s~qN8BKb- zt0_HdjLHzRVZ}o{xvf{PSoQMCeJ<*k+y$|O_PUlaq3b&rw;vRFk6(fn@bzAG_qd*u z86ito@ACtTrX~HI*sjZkrxW8^mlI1(u4cqTykDO3xeZnzRIJE3kj~DwsXUy7@xxTN!zipu~^Q zD!1$YNC~uRU;Z+qpy(9AEMm7^t02<1OI{U0;C`p4gX)#p?#J9moyz=!FN9l#Ft#Is zY5!)qkQ4cap=I39ZPVRJ#ef`HDcbp~e9J<>?`*~d7cN1|&T&?g~@JuwbEF za0>OGxbk{!0lU@adz6HDTV|ulcdAp+uE*WeZU)*G9va_dbacVLKjS`+%Ph zZa;W#J9-N$JZK-PL3KegTIMih;s&%Za?QEgfb>a%=$&MxEFyE1v>HPV*}wPa-~#Rbr7r4)*T8wN`T_ zlM%ny0%_XEN|s8o)~W90kaPCsrcawFrHSA>KneZHecok~BOnMmy?Ec7SQ0M_5t!AQ z4srFoqJuFY8c|$ZO%?CcE9D3%puLEw#LcK{+wP~37dka&DFQiL!r0FUH680Sp`hE4 z?u^QE*X~tBE=0=j49GN@Z{ZnFLoui!kn(fFMx=~%i>kz+ZmfgPC#_N*(BK+6Hqv6_=|N> z#@xGmk&hy~y==MvT<8Zh^gMNb>acoR?O*+2&F!?9zY#@;)=~OvTPzd)4E3vLy{7dJ z`)^DHpQaoo1PDrlfPw{>&F>#{(0iBBVQ4z%P21y_JF1EpT)@Fy#;nLh0^fITDO~M3U?`Mw*E$Mw&MA@)l3bj zQV$Qui-am5W=0M}2ZQ)*y@m^Y|a4 z{4g)nkq%lNuE`}=fC<^I(y#WRYug+4Kas;^C!oMv*kF(3rA-N_F=xb7(`FYUpZih~ z@nPm4B-3z__+yC~luT%Y0NDLA^D=Ap)qUF!ucfjfX9qG(tHWF=Cb76sV4xH;9tL@^ z@WIlP4%IsudN3aK_uo~@N5343)cwivAo+~SaevV3C7Cuqy4CxObh;7mpV(!u+)0Ov z^>dD$s#N;^rXp28Y#ZdXjMx$xAlmzYwU)zK&2Sxiw|*k*hl6Aru)OjxasJ=^kQ{Vcm+`<#ReEenTkt#_iuWmdm6 zp&iCE$Xwy8w*xKDo{epd1pE1QuTj%P@46MPXDM^3`37rL!1#Nd+Tvu_+%f7$NZr9C zH#Vf7G0}JatZt>l$&#MCmRny|X%#VI08-dR<}x_HDw;hJy;Rlr+b7-BQ^;%G_)7S^ zF;4&EOoU8mech}A+)LwRVLL-KX~6@Y%dxXCQ5L!;r$KvQnOJmsay!?%v%P%Xz0Ou( zQy{Ht>6-Vdlz_XX;FmqoRHsbai64(x6ALdse+Q#GYadTY?#xEN2i65*?d)c~M&4$K z73a*6+=Ev^uZ}*%q#qEINri-sh=~9LFXA;lCyJq`DLU|DBOp*%cQyfxlv%mHG!Y zPk4FFzW45I>QO|DmpR0|v#rL;BA87KD7HyImYcWtvlvWdH`c{N_MUGmXlqH1&O6b2 zdlF;Mlp8#ZFqa{MfQ|ekfP3ktw*zLx-!`@dRVES^BAv8+D7`0F_UL)ACO(j-;`HUS zb?Jtk>}0?)J>>Y}|BJ15-0l@dBZMcj%ytG<;YYd7TC(u3{b7N9PK)l3f{L$o8cuK} z`=MaFI~&fPcjvP>gT8t4tQFOgOrBp8@Imrv%Kxs*9y&3Cbl1Zn+L)hM(r)l1$4g_` z(8ndTvGeX|hS4M7b(d%`l$jr6eDLudAw@_;2Ti^H=cFsC&zxi(-^*0;L|xQnN{~mO zG!vm|)xvSc-p+?D4pbL2h3=cUIZo_3m^lK7J$Ua**SyH?7R5-OZIPWZA%$U+FXgB+ zAw}=g>$O)rDnFCDHG!`!6<=UdTnwcn#bzn0ou6vzMz}REv-C_%)-ac|)mh&18GGv}uFH z;O;?=3S*CF!ZiDW!tjH)9Qey|1x6KnDv^mI5hg>j9{hrf9 zXuOtHa*LPSX<2PzM0I2*^sjJMTkMFJQWg_~EOV-umU1StmH`&(RQ-@p=3}abC zyLZOyJU^%$D;cFxvh>8cL8aYeDN@_}wd!oMYHeSq_&$AfgW|c$&&5A>{-4AL=w1z< zq3)^CEgH%MIu%_IHNJ{znXNJuDsxb14p`mr^;CwQn&JkO-?{?KVtpbPM$0y*%YY6y zTV*#}RVafi0K51CljpW&83_=psY;f%xyW2L&{e+e?881{sVOVT6YIzn=GR}K&Cy-#(VDuZzOa8|XD9$N@ERD4_BX9e;RhmX8%#S%(M!F;eH>H1i_JR#pZ_W7f<)ruKw2wg(6#?W38}_6iEdwedv7hMyf1rAx!F$wmjm8{U z@zTqgrbG_YHCkCN0S9M9)QhyIHr7W3$&wwvUiu3g?uu@UI0*p&EQohgmmrM=^)hEx6jFw=dWhx(z+i|(H2$CT$>c5D}hZ-{V z-Zy6sB!MSu;!X#7-XUd@?>f(a%^J>1t=<}-QMJ5OWEu9KEh zjyo5VP^!UFO$DBgY6xKe@wLjKBV&SHQ_GdPQ1iycGe^WFKvO>g>oh}~`L-P6Ef3J5 zBz%I7gxgA$U2*%{q$I)Ln91KG!ZDo1 zA?2^f%NT-swmd}sg(9y zh21M;Mnz;BK!Kg)yW2hg;pdN<->#laEkNT`9pZxiK7Tx%)9Cow;z;7NYrXN(Zo5c-PU^EyPrFZWA0{OumIVYQbK{BYA*#C6%?W`_vMZ42-oh`9tK+SqYctFI}viEGG{hc5! zLT+7U$t-=ZQ*VE8)Fc!iuP7h_2-wSLpH~ zrFa56A7X`Q@-V;YEe&M{Ha8`Qru5abte4Bk0_KPm4otE9W%puG98JHv7v2cwaFoqM zdz3zZ2w}XUN!rqC>NPqJQy&oo2?lMF;$ll?`fMxVgBTwa>A|7%Z1 zs%uFC)0i~}`#;EPHF_eZZ0?%jI zz(RJ3^DV>T!M7?G9?+Tn&RLkrOUfm#2VQwqqP@h&&~udLv}N9FeA-um^_p=}Ej)j> zM3cax%oY>-aGbr(VFGxP%szk3r4Ihi?D7idB*k4rl(#5B`^e&O)y$rs5*0nxr-xnF zcD#BS_n>;|<=W=6mh@bm6Mzx7rZaayN7I5-A*~$o-(k>t^UD6j$`p`$rNJ*zEZJ+V zRSf9ok0#vbr+}0Q33*7-{M$Q(cu464Z3W@)ce2O|xQoM`1m=Zf8WR3UmQ1t?y> zlYb4S7hY}>e6XC%z*$+gxi7W&?c)>6RNXWgS1Tx<^~mrZ59FIZ^{5Spi&`tKA5$w| zie~%{r|lb>JqcJquAYd{2fU)YNJIwg^~^4zoBu<#^>BEB>=}Ibiw?7Cw)uW_w_hhv z@JSvR(V`3D58b<-kKx>T{8!h-$NgiLY>G7((4~>G0`~lSphqIsgsPwVKBmAXvX8&w zI|1UlrZ7T85QDp9b2?R(gDs5lbY+4-mMW_w~p>k;c&<1F?K2LylEBh5DOM_=5bo`ABw0VNvlLrQs_mleq zq|!bKmSI9~H!f_GKJ#RU$cBH{&eAIxifStNy z2%}mO$^-AU%Sbt_a_2Q!dnKm08(dyaqVZ2 z{N;5CorB<>IWY}*cZT~LC~pDxls&-qePDE-eTlBG92!dfYCi=`N4I3_n`CDAwa7{X zB}bq`5gG5^zNdpRaC!-P;$?f^7vO+ZmG7k8`EZ|$Z%E;T8PFq-WD|&j2|NVgT&IxD zg;XR(9N6a;YxVVgMW5v#V8GY7KmjO;2r*5 z_NRY(*XmkmR$rI9{y-lQWD$h;LbY#Yhc`>os9f{Y4Pkeg7d~YMKI)rXo#8VZ&PKn4fXoGDOOc-&fZ8aL`V8XM z#IIbOoQVhLQns4jZ8~-0JjiWu8{Vam60qa%+qG}iEj^{C)Eisg6k~w?Xku=^pMs(4 zDY|C5*lfpcb-+wAVi+mV6wQQ(sR>4OKQ}gN7!1uutYr1^kRlZIUm*VWENiZS4nw;8 z$Jko@cMO+rSD+^ zHBe81M*YEbEMv$|E$(g7Wp{sF?%402o?Z5=#V3p9KRY1fD4y)Ghyyx@xg#SG*NuS> zcQdVkXWS-tcU;Wwfn>0P?5wH>X5_|1znc8u+HnUNPlw`Lo8=#9U~gW-%)1DT@-%wx z@AVb)f;xprx%zVGvTg2oUF7ow=|Dybxv>uItrbBSqYeBC&FcNJF7e#rX(@_*q76aP z`7uGby_A`p#HgqNf&d06l#vBD1L(N=foEBigfo`}LadN);4JB5i6UOOwJrOi1A*iF z8oOI@2^K+HNq;`sDcn@#(b}2JLJM(w5qiz_>=V=G{F^K$aQi1d(4K#HpnoqX$U%KZ zOv|9avT9j*oTo|fKlMnlJr@_(?XByWk7ANOhIB$7Xv}G^q~wE00h=hpgdJ5CZVL&{ zTD}kIoNX7&1&eP=E{Q(Ld*W#UtbR*X-oy-ovk>;&;2N>sIf(J|sh(^V5q@90OzR+6 zRTx@$+47XWsOD+tuAa4RhRF*W+sUV+lG2A*pT03I!lXC}fRwOr!y&E=mH_(Vfkz(i zfk$oeWGQ@F(H!hFa48jV`ipPa5_5)FWX|*9_p`hMmfS_BP?kOs#P>nEs@<$+#Mau@ zbiou8kcSKvB@J??&NphQT7h`4E7oDB#*W}}6EyM>Z3a$KYt5``OHGSDN@986ehbYgQvUdW0X{~~vE3NOQ zfp0bqrOEqS%`g|;fz|GLfgRu|CZY=~%GFl`V1_Fbq2k&ZRZ8A3B$Txi_Xn*+YrtPX z^S~4~cHc(AmI1s1I#_-c_H{m9$A@R@lalBzFng4YE-9z(YR+Td5$*LSbuAcsj@xi` zKyO?JdzxoOUfuJGQL0uhOTwr=yORec|BYMLwJRX8EJiL(U)4s-aO=gkSrX z|21}gjJ5FnaBqhXvkRokJu1+}RZW+pU47iO`}1sy3Wo?;Q(+LJ{yym5Ez~n|U?55D z!pede-2M*m?m)vP79LHcG;8rm!NX#|VXG}@G@- z=Y5HeEMd4;X+GE{u><8wFyt>BzD1_BRoGn(y$GPGuglK#stagVNs{mWb#`i>8)Q0FXKz4yF?*l&B&0w@ zOk0rLfT*fu&f|Ah(#+WM&fA-BCHYRTijTJKNQF#B$=MVUF+{$kxtg<67h%jF>?I|n zBZ(cAOj#;t82|M`?)^0NrHysH@Izw?F<3Ln6eoydsVjR4+kUWf;FX>A441^+1pUTn zYdV*#WFsViyOQ@1h`|vpN)$77!}14%ZGAb~%I|Oe+Bt4wJDn#>ZO=tQfaaZP$|7%R zH+7-51&_X@32?~`)RjuUNz0`H6*DS2HSKSH{Inn#GphhYYS`_f?x*yY8KefX&V9bs zNxEsQeme*{1?Ub|)z!hfAiRO|@T$VqP(~-Gr4Wcg;@Ex4{lTiYtSAZX=f=jNnjfN? zbD&Z%=$o9mTKLf2<6$Bv0IzQxPIcJX70(k6(;!0;7?s$)t&W-h;N|N)nme{$)4$RU znH9^AN3V*@Fd6*osE=i#ku|^|aY}TGFtGfeF9gd0fRb(Cf$lDVQiB{{x7|th0 z?hn{r)&5~ny%oYr@(1!hI%OsoecdwfT6?*}msWJ#vs^>bDbd1El>05~#8<4k+>d&HrLeo&jQUgR}tKSBd2E2d3p>^Vn)zgpAKf>)fH1x{a$ zI|3^(xjp614A6)_EA8;QxV7gTCTD!*UZru+Qk5!C#oL|ci#99ofb;^!$3$cUjxLG2 z)wYH=tvK8Zqn-G9zDzT2BV?Bc#LvLSdRsXW$Ts`Ea&mT{nRq`it z#TQ%J{hOna|8-hEP>WzxRp?`<|7H0?viEsi>Xw77_*8SjnT;gA^JhceCXvU%W4Sd-b#{S7%##(LB}vZq1rYRN4NZY2*&aO+|(uSOc3 z59l2~sn+TLHhu#?emj5OnX3HmH=1c0u#b&0c@fQ@?JMTi4x`SGqQ` z!waMx`nh#+Tfi9m%iPYuC|!5&6b1)UdNb^6+(A7HcIX?FKDRq74eECXG=2Np0V}pl zeR=c5^)sSKiKIn9jO`129p=fet-8UP1f&3A%N~JJp*IFUax=%tqZ1p1!rrQj{!TLi z7YmX)E&K8Aa2Co<5)PJ6UaT1IkN~r&w)7OcA96L?9C=wR^I&mMw_tB|^ZB}JviH_7 z$&8(V*4QUi{};=q+%A*BY6i6X$M>Ovq35;=Dd&%+cZwLy!Ap zIzc?l%@4^D;GAoU_rZ0gh2>u(i%4t-#!B#Vd>h9Q$Qd+a7h$mPKD4^T;Me`BUGb(2 zvA?M9AiH8`_CU%ek!Tl0b_XiAGfjpjSc51TIfugqGf9n%Xh`XS(mM^`FMza-fcAZV zVfOE?Bd_P#ps8EQT}sK(=5;nJHkZi=H$9}+8=Jf`Qk@Fp(r`4dqt?yARh6&tI3DV! zpk8im#Wh+2IXp3n$S%UT(xt<=&FcESTQF@6QZCyBhj@b0Z280zHP=gab5x%@3W9g; z@yOmq^bhpu-rg>!FjYJ3FIkRE$_9-YKrB3^r%vM8doSJA6OJI=d0XQBGDmOe_{H%X zRVF`-AJQJL`!^>r15P`{|3N>A+zp1aj4`?E5z z(D^loAA77M@O0(}cG{QR{yy%&@WiTtn)L62w(s+3Er}E1GtMrZ#~h|x;jN>$oUm<; zDD>fC(>siC0T^s|~1ros@8X?-6sJ2A6UtBt#b)va+2s>9VEKT(1o~ z{;fM?D3%UyyhD=2Mze!CvkY+RIM`chL9#6b|04^ZRm80aHwmPW1^W_pG)nouejURV zn=dt&&?71ZxsPzxIU^9>WDL*&S6x z2}O*LN8NJ7Q{?uK2^}kq9T}ShZ6rFu?ne;GG9PSRkXdgRCYB&*6(gNFCu=AdTw=Y8 z#>qUCS*L@`9BJ`yRlD+mZWAm+k_rnT4Fke6-Zyab{{yZ@5}SL3cv><$_iI+e?E7B? za~D|abD1Yo8PpHnh_h(asd_@*1SKwVg(<+PcSB|H%&m6GL~}h%LVSC1vq}d0{jHc@ zUF1ECLFc?(NA7Z|Ed#QYLD-8ny&(kFqDTMBvDn@IBlawV1i5!0Su8OcMMNR?9Sl~VwG;+rfwqU}H9O3kCOra0dHZ=CJPvC-rGMV4KJ#@lXB zLJ~H)@D5Z37$p(1Y%sg<4tyto=h|_tI>|cW_ufJ1n8+KDdWs)eHle*CNCuskBV!5V z{DQod4Qi&(s;+mwW#DM`#yk8&dH62hJeK`C>)DeomAgKu?D+$h$upZu+~2|Sx~FIcfcZZ~9}poUF*DvBI0OwK zvNCtJU40}IKi=G?E7s_<>ELoW7e3U-+5oA>^pFk{CLg!n++B->ZDyR!!#?2#KbGIHrw*? zAf#+T(@A?%V5kp4*_D2aL?l|af850I4%&F$x}3Oly471_5YRz0At%LINyG>#99SE^ z&Y~+)BcKUMW4BnAY%A7zzE@A$wGL*=0H*iuDpDTDt+$^A%@bm#f(y`wDpT0v>55^_ zd2RYy`M4jT97GPwp%9?e()xQ|2Qxp%)J%~h7R!dJFw_eHM~vsze+WIx_|?*6vAD+{ znv@;iX7Sj;G)fmoRMb_2Aro`%dKN4&0Ie;@UoAUczgwWw-;957!~*dAV1VqiudoKi!cJ z#**1#)m-u@RFk>x^>Iv(1gF`f8VXqh{z*$#yW0^lim4r=?r`tbFx&)PA97VgR0HT% zM8I>fNt87mzg_LRTmmqA0O@WBU66T$10#4l%qk6sD$~8Jn*un!=<9TaArQ|f#Am4( zp&5qIXAo}w0XsKlBmoZ~r>l8tf3Q2(g-<2k)B{-}FQgYs#b|8fNsm5`o)eofvC?Or5*qYc_DSiqL-zm|^O{Hz zjGA3)z0{#;E|uy|MI@+uD)r6B$Lz;bcrWDB?N%hTo}=hBxwfo(ImJ4LkI{om)q~?_ z0RMVS_ZW;)2Md-Ang}V?Y1P^Z5$EX&2cQzk98e*Jy)O{UY%j6S(k?FRa}x8k*Pm z9_kSd&@|-EVdF0>k6U{tpKZ>xBZIP79?cab8L)DJivpW6HAsPy0`@kS|xS5o{P~}>3oSpew4G~O@G!RX)-CG3>y1;c74G*7@rt~IO6=F z$}c0^{{^U-s)=-FolW-PJF=7n_a4bQu=dE#*#?Qqi6)cSLr3?akTU~|h1A>|epjjR zD#r?Un25j3`^H-er5&31ecX1dSuI_N!*icZ1r#dGKR*sV%U z`ZU(v8jrQ@nk1x7dt&Cg-to)Np8>%<>qe7vm3WCd z0k*jbdnydrRC?8eOwi&nsfA)>WJ886adhKq;6eF6FNV*BH%|>Lt{Ca4RE|VnKYo_y z478;{OeCa|&U&Q6%8{8mb>|%~57r{`6J<qcNclRH2X0JomXRszQ1~HS5*2q7v#b zp~@?pin)r4`$C~6e6(4yfkf2(cl@1;UTGgi1qOhsK(BVAE~NKMuY%KGQ&nn)`^|cGtZA4@>29Q0x1I(V? zLL5&y z>z%~voYuppKW^^y{m1p!d`9h3m2FR$CwgGOWD`7TVWy%8@BVS9H$je%hll(zfvoJO zKub9~WGV7roz_bvMW&$7eHQ(rh$htQSp+w@W1|*446E&!K8keY@Sg$#KB$s!9(;b1 zn?u8L0kS@7oCl>?dw3=FqdOF6ui{75P3@{bG(3l#39?T?@9WfKDo{y6*Qyb`ts?dR zalC@q2n`1K4~`*pQI@JegPg4qbzLp0##I0*&*V^o6JmK48E4@ws)O_kw4#d=uPHld zg3R^Q39+%Kk=>=rk9fa8Xw4@8_cfuU0xaRcl;?qz)wODx0JNP&NYERkCJ)j7vl(yt z1XNdi7Q2%Vq+%>+gZ_Beb3A(8m3AYj!njILEWB@ZK(m}QmDF@7mB{X%sSyx#) zzu+@3)~OTphJ!6p`&AR9UF1I=rJFB^1}oe=bQZauE6LD!=#ZW)B*qbPoy8Vm=BM@2+(3-IRm-W>)(=Iracitg=;O z%a}n(l-TL{EJ)Me&Cw()r3fge=+;Z_Uw{{3JX*a?RRRJzsf{dp+N+F~z?xV6e6dNN zKIo;$z?L@)1afF>DJAz)JLrJM6y7Vy+nlpWX6{Gzj>}W3ec8*_NPN-ua8o`m?5}E3 zNfj&oS!lwTa8kLke4St41u-Zg;Kx8-NCTr2PPWAJ?Qn^#A4hGduwXk!vA+&N5HJ&C z+>cngUYB$^6`F$_w21Lk@#OujdN&CimM%wC7uIwG>s)20WOfdOcdG)rGCV*Ew~5Gv&)!hpP(yj3VU- zk2w1B6Cwzmci&|P;#ORPe+lnm{H)pQO6+UNln((%pMVJTq)M4DH!x|hUNtl~2ZHgg z%Oenrg9;f9R*h9V2Fa|ddIngwu|GUHT?oiArgj@_kZrztJAxY@xH*C!B!2@k^Z?>`F5|N!lud%I$K*1h(T=te8yms2TKxoxDfg}wh7?%C zLBX;0*&06z2CS}*F%hm7_Y=LevW89}SSLrvsr3%X`%~$f;`IH*Sf5l_1QK!3+A#7R zogjD)*E5bBEU(;{0Oc+;%jR!KUiIfyFye4l_dgjVBvmldA-{inkf&ha*b(fFRy!o` z88#6lrhElVRm;UNiUqt&zyeCDOPzqEHY$$dl#Wid2Cnf~>E5w0ivPA4j6KCmG!mp; zee??8sS{Fx=HFs27BhibPl^JNzqGNq5jONrjQKIOcnxK?a9Q?0m+x4;`XSJ<&uJWR3?T3!%R-_%JR0 zf1WxUD+S3>NJ-@031F8a!qdbLfbZ`rx|y$fY+_vsri6ijuYF030{TZM#PZG_nnTtP zwq3t61<8Q%Xfv8RSQqWr{;u5O+ujQk?xXDvd-qNv=hwMnQz;$8v&Kz=x;%@JBVYQ~ z9S&09o`zCVZm_VkDm^*m0>`+o#1x)dt>yQRf{eQd|GS+Y=-=$dcj(R>-{9@-xceu)L#NFzIp3L zak0F&Q+&JsXR$I2A*Zk^HCl6{o%euD2PYLa`+v3e<>65F@7pSot)h)k*6eFj*&<6A z3_@fWQnnF_jN-9{?8z=$V_&k>*g~jBWtj>?gyea$&m_-deT?wCZhd>c?;r1b9KS!_ z4u5nw#yvCl{rP;Z>pHLVJg>W;fb?VW3*PbThLgK&V5602(JJ~0oP0Z~nC3${^fNq1 z&%z_>_r~<=6!e8mC}EBSo%yp+x&yRfs*e01eH%GC=ld1hphUG-?=#vQTo{5#ByBTp z9|o%XpG)yQ5Q<=*a9uG4(Ko5Z!lA$u!^0#<dsbs^}Z=|3RV)@rnd34vp=&e zc*nj9IIg;+WNyClk@{WHy$BqM4V%PgU+tM=Qix)GP7kUR)ofc92v0!^TsyU>upIX5 z(4mfAJcbMB;600KA#q7MCk`@)fp+=b-AN8VfO>Q;4Ghho-)Km6gqfC@i1qe%$)}(H z@0-@D0d9=}vA`ilS26bX_*eJr^p8aulii2jI3zIj`CohV?~`Z^r1ulwBJNhjL#3#{ z(*rp01DAN4L~+LlKC|Nsn2!m zLN@DT)*56Ic93Vp7n0S1Ynq`vi`xc8p&&#>fuO$YAd$q*5>Sdh{ormmfOckUGXj(O z;xWrwNF;+QEEqEACwEPZjrThaFu?*hdvy%RRArSMl)*R9wQ@qtFAP*_#kR~8(hUgD zt{ADCY{!8O!GL5QKVDg2;lKoADwRG9^+S;vD0tDm;60C+3kFv0fO?pzM%$_~LuenN zZR~lj4MZh-b&o5miKlEP)onBX8&dN*!6gT>Fhg?mrMF0Z{PV8-_e_28w8Hji0SKuf zUotk`FkPi=BhC@7DPX$^-xf{}ZYEHu^?J*)A4hH~Y8Cv=CzNjC&UJmUl}%@+U*~(< z?#`i_z{&WkR|CmfTu<%fH`hT~xQpS*M@i6^R9f=h)eGD+A2Pytn6t{z%JBcJ4Y zq%x#xt`&2EJduel!(Q2W4@_i4b{soUR?PwU0-Kdc(5gAOxf$4WS1B~tOh1!mLoEMyF#x2j3sL9wYer;a{{~~(8GDwH z(dB0BhAVjr@>FCXJbDOlbFWu>kNwdE^Uh%YHWN?*)`Y&B zdBg9J)Cv0cKnK7N$4}J0h(x~gEIjgj1^8ZwZedpxgAPalAR~y$7Y6|tXGkaNy?9K% zqdEk4B)W=dL?X5hCj!?&-0;!uj9)bv)VVw=<{7Tfc_&5Vlpx|{(+kA6j35gdy!_<` zx2^@WTzmS4xOp{TmosmoKp*z*N}s!&&VO26_R9|2nUqt73ssQX!342BK;Ud zmeyN|J`Z=S!$kae3Z1}L286kUb!!&>@{f`TQK>;CAaz11gtWgHf0mHf3ubQ4W^mmS zFt`SLE-3!;jrvgzqrQv7RC+lLnIeHYglzA8WppSj*Dcd#>wVB4*k;QxiTLy_VvEOA zT});T;(Vc9iotHSsPTd-$Ruu1oBmYV4+J;Rxw2lhD8CG{h`dwHNTf7qN6yao%5I_> zmELl%7|;J%h)p{xKC9RK^8Rtt*`GchaLPev@j|wUdFcKf)o};!B+yu+#y)VlNO)6!{wS0x{kz*(^6p7ws>~x(J_OI+-ETrzW8= z9m3fn3W`x_bKFVm$0HAPk_@^RVHrg^cWVt1a4iY2qWE@-m))rc#s`=?V1)5WeSL0f%vhD zmkN&|1rpr-AEMfNhGY%cQfc?Xe%kCA(k`%+l1+!)ZUrU4O;jgFdW>QY2s3qoL{9Wt zqwRP9*QHQTm47@aDI=Fx8}ZFoxgg|M8Y2#1WRh?vVi6AJLd+xQO?Go9eS3aIIv+X= zkbJ{``dnX(tg0Z5v2j_u}z9A%1p0-s;jZU z!w4D~voKf1OAv1ES#xN zKHAvQcTB(jAitik1H3A^eO8{~*qFr|k1{~;O;|kJ%12-6rxp%P$3STLc zIRIj;IZCAnyEY;Rw`5a%7~p(Ph7VRPj32nM9pT00VJ$|`K)mB_h*;fMY3Wzbmi_mg z;zuq&ioD3d)%lcS%B5)Vftlwl7JGFEV1SqhktrcgcW+S^g48%sO$T7w3xmFHY?5EN zkKYBg<$5H|$tmp;{0-(&z&k?N1~3W`Cuy?5YFMJ;K*RAk*M zbModZk%*pjPnLx^ARO2`9xg~GoKMz9QVuvpp3#A#nqKwNZ+`1@ZV(ZR34?H%_jIq zM{wN9%~?H6HBi^F#7Z&Bfc^ZmfM$3tYfWbRD|vRL!(eic{n@Zoi_Lw5!RwHR%^pnY z6rwQ|m1{S^1NW9+&d&-|`Zr^|1QTHfUT|=Qa~IaKRzs5yP>Lz(5H1E&z3UhFvJ+w& zV44J~AC;x+b|I_hy<`J0bT`_OwBg2+zcyOOYvU@W<{T&4c$D*b9So%I`x0e>BfF!x z=U4v}_!mpLV9WDtVT;KG20&e2rV=1ePDD&J%^g;|>ba+^mlR zx4>{H9OyE+j43 z`FC9D7hnUtmI? z5}3C+a8uvNW6b3;`ymz-neAc#L*_tl7`p-i0P*c=m!BR}&tjW=&a4@N`1XCb86*^) zX5SOxF#G1T+m(%ws3@i|86>fEM>~N*7f2bwFp730H{;Y1S;T&>3|t9-Q~+@lJi{Eg zA$oSsEa5II;(cH@jju5Wecpkj;rBa%Rj#oW!1bNO49Zu!jf<_^j5;Slj$Iq7*Bxx+ z55LFig=dA2XP>dhx;BBv~YYD7bLzqF6QXRW$0w)4-?EV^9eSs>42D z(D~?gnRpm1?>gn<2JmA~J#6HTa}N;;!2kJJk0CwiEao&~2j9LEArN1Q&rHN;f?Xvh%z+!mKEY4d&2Iu#>!`)i6kmR@r~P@^N8i`jFejl$3HT)M zeM^q^yOhtBp7k2++E>!{j@Fwun*dnXE>HK5ISg5R{dsn6^;hvl^E00Dfz>f@^6dt$ z9UQxL*F9P{8CK*^iZGj3Y`!KNf-$^8y&v$r1c`^OY(u~rK^u$j{k8)SOiQ0Z2IUyu zmtk)&%qaaa4(7FU5HnlkD&N7$U97tuhoZQnJnXDZ2!SNLGcK>v`hf!9<{+XO7qZdLMRVQyJlD0{!zcpWHo677~*{$bP@X{Vq}%MKj0K%r)@oH1DNWQ>Z?ht~XFB;OL%x~sbv@JB zJe?IEDuKD+F={-W=*dh~TMS;e$t*Rj@v4)+x7a;|OHeZf#yE}o)n z#gJYat;v&OWIJf+aa3{v+RqqTf0f~Vcbn$vBL5o8mU7*ZtptyA?7}d*s_&zfeGcl? zP}W4Z8cDh(qryEMH8->t{bGnPR_wRDS&%$o4@zbzEC#g)?d#E7)!#?m>PCFxYd0#O zq|Iht3|=`;cPY%t!Y^-1&pi4z4hZbGW7@Ri!n8ck$Zrol3jujA>GJQ|=CqIXK2ef$ z&bq^LwbA|hxeeQ@BLw?yI%08|eS_Yl!!wg&fD;?;m>-OGtx8;|&`i2%J zj@tQi6Ra(_R3~@ot3IZlEjwbLA!g>OfBzfNcf;Dg!QbHOMfK{TB{2^T8ab)QeCU#W zN$rgS{szoT!R{-O`lAbrcVKq%g1S|@@WurW)>(4zT2br6ft7(FiJ0_~lWh5oZdX#h zH_tskMqP*us3sF^tU~FFs)gGToeL2GZ7VhV2l%0l0TmF&W49Z(uA8~rT#YBUu42Wc zJpIX0DC-@hR-A_^ju!SV;!6gH`fB|~=A2ImlD|)##~by#W79lF_cil*KjZH@BZ%XJ zB5af(1H2|oI7hoHmf5eGk9U2PFjmOZ3L~37x_M+BHZ%|1`~>XF0S=*gCBr*uLEHxB zE}H>~-&r5%7SQHr^+r@6COfu8I(#CH_b)f~+#i`7k&i3vX$Y``sYmfQN0dv6=2U_t z&eSmn^qbAm=lYabzE;^%a`f9gi^11#7XzfgE9el))H_ner07@<>?isM1unb! zla4=#P-ih!x8Ge)Z2}Sf7?fXr=3NLewwA$uwHK_jh zyYE9@Lff6&tqr5W)*CMlDSt65>%^wl-ZtQ{EW8{=$BE^8U3*)g(@98P%@K!zEig=E zHEjq+iRAW-6pkrA_(w40hrNERGLXH>Ed1#~GDl-E!na=Inf$i`bHDH8G)s2aH<~Bk zYgXPWd_T^?L%*RixjFCIOL^qA7U6s3_W0WA$oG7QSEe}wzZ(uS`BB_qC5j#~I2Rf8 zIDz+Xb1wpayZpOvp2DN41b6>N!4rq^2H%v7$-cF!$0*C*6^?Cd$6r8bDZ<}w4t#v^ zJkyT3cs(u2*VxZ%O|)G;=|Z?9P3h44i;Mi3M+3~y|M;t2rdrw?Lmjd_7y7(8vR0is zy>VH+gxFKCwFApvmz#gv){jXvXjxB)iBx;vCj}D7Q>iv9fx4eLgpGKh9%Mz#oZSmqM zvN}L~4qQ&E1DkTRJf>YD^0~}W;V;_29<;@3(>8BG)i2c9`sjdcO}=uJE0(rc_{?uzFe^<&6G^;N@pM{o#N6 zN#5JD&Zq3Ma{Ge4hTlWm8;irxb0EDs#U6p(T83|Vr!gVn9 z##2!2y8}5faBpY7dU^j%c}&>ji-4oK+Lme$pwqV9)Ru`n2lT`*E>0$ z=c3vTt>0C;x;@^u@9W;KQS>2`K2DM8WS1`AvQgPlhYs!Q&Vdo%4h-+Nj?d)8O(!ZU z$CJMKZm5~*S|{`*1Y#>K3%u2A^IIQ&T%*!JP-}1fIFJ1m55C=|dn;jm#g9PUidIY$ zb)s*5C#dG%wHfH?P3I|Mk1aDX*Dy&d3aYP2ZW^OD<;UjzGwGiYSP-l`@TRHrtYeNj zwXVE6+4d0qx2S2ZMpmbTmd~rCD%Kr~*xbVpXXRZTK0jgciWS1I=b9v28aiaVyfM`+*f(+1a7T|0bP5r4 z1+58cANh&C5RtYx&IV`JuRmnlKb_tE^FTA?P;lcSp;+_)q%v9WGY;Ux-A}BI|G0&Y_iVH`pa1rBI4JBGcf$w#e`WDM0oATvYsvHL;CKt4n<3=dvrnk=^PBi(r})!6 zK}P)T0ZopN6nJO>h$|mYdka5LmB*(I-SZ&Lb5YINrYkP5ugq!jQ)C~s@+F?xc88oC z`g<0V4+i`RLbUgS_-NO|PwxY6_n>a4y#xinWNM5$4;0c9?#!ADfb0 zod9PgKy;v?5)8%Ma;q2kMv_b*omT7VcbT(3$OzC?70uaFWJU;tEgL~rqkN3JfP(;S zDMi^p$kWw{3(xJm#{e*USKbPHzED)IA6w+t)j_fD^wHniHBWz_qIP6no$=O=(goOf4e~>iPI< z)J~L|51Evb$^IwlY8oSZIVPVw?HvIh)DN*Jux!ZIdjib@XcdCnntPl zav7&)5)HbW$I`8d!c!nFGP0UUY+LC`o2Sei^WG}Lgslb%!Cj<=9;@t9L@Fr6m z00d%zUto~ygc4NYWiQK2Ha`Xpxb1~h8CXb#TB4gn9PJ?h8Z1&eU3~qEy0FEwC|v0% zfOnbo7rXg7(DVr!?a=GMcMp^z0ZlZKyh3R0idiRayq9MOmdewImaL176pCi&y1n?}q$ML<{jnL_ z#M_1%Ym?jv)ijIV%IW=OcB#`ZV}(a+^`5_x)A(%1aO|GF_qFnzGQ770)VC*X4Qxt4 z8x0Z$Kriq*IyJ@%*TGL5KLlib%Ud}(6D%J0E`L~QX% z8$aUK2|_dUSX54U|M=MUaV6EmLc>a34;TM!&3DEK*#;-x%;=T;U{S=?TaDeJD)g>J zN=iKp1uOF<#phL70gnSL2kG5gX`4=?zTPslQT{ihRkAkYCGduoSwY`E`qeDzl~O*O zcOAMlSkJQgNLOO~$b4E^bwf}KV5tMw3m+H2SB*wW_$Yq?lcgz^XKMFkfKQsC*7Z#% z(}aA;z|VVxcIeMb`RbdGv$Gjf1`knQ*pJ)gpytxIg*}xMFQ(;94=1dZIV}cTpG3GM z`%}DOP$zR~QxFqm-LQ_!xSvmws(p{C_` zrqsO2kOTtjS5I#1(L--#4mqyte%8{&K6O;U9Gdc0Zj7)U^1v&E-yJWfJ+7Clc;0-V zHc)De(^12~ka=_SN$ou9t0mK7t0w#TGC6bHO?+i;b%QMoMym6<*|!b!E6fgL$Pr=@ z4ghG8*&MZ418M%gDa;2ZAe4)MNv~MnyE|<8IWTuj3WJ?Oe|%3jw&K3iVcO6DpvL7t zR4ZqE3T43b0k|V_MOon%)MW$wh&Pj_z31S1;?ci4^Y7dM?p>MIXguRX$8kuf>c=ER zaz)0D3>efJ848ce?1}*;KtIYRJ{Ag-lqxBevL7WlW-cYlnp5?zFu;BR^7>Pa8P0!E za*C{==?B#nzJ!o;zIIa@1(u>4a_k5lXPC2>Jrd?wY9iSzcx2hQ8f7NYr(v?l`p%;0 zV)cvIibV(3iLGUj&7m!OMLotS-AwQ0<0k)EEkBoq%0b?VbU(bff1fqVx1i22SRIUg1ToEQ%a z&RAi~!==3+cPnn*O^NWO3CyH2s*3dH8xO4&1^cAB`cX-Auxg9{<`?qU4RwGe%`cGnrHoIg9-15(Y`_I~ z-N~3x0+!?H`k>7nSOtOz@0O>hxW(9qYAc3Z3!{%>zx&_U3B(GeiZ$L#D9M_ayxcFC z3Il+xP3o*0E40sTXyPoJ*FVmm`MO25Q`>zXB5|ecz-EDg)FFI7rJGVZnrmh3R0y`5 ze~Jqy&sOTFpl#B9$qDO4q;cAMp8WyRc+KwJcp3;-iMlZRZI8LJn!nERy1G!hwKb7f z@`&58WtlBm!JxTyAg0k`&?%?zW7v_Iy9t3Abta_JQGtzuF~v9~lkgXln~S0KUr0tP z;{qRC1Znerb|e$gs}f|>wMhRHo$H1J)eaVJQz~iSwH9p-zkYQ|B6%TK6AT;X&y-7b ztUEb}ayHb5QHqSQvQIXksUvc)Cu3ibN!EXIZKMbuVCZs}<5?5k)0Wb}J9VLSK>@ld z=b^`{$S!PMto}9pg^A({&HAsY7OFRH2-Cf?9;vCT4xoi1rp76Jw{^NPz#@i7)f^W0 z2@1a~UA~0CqVL+S!m^#}$$?j@hq6}7TdOXqVebXRze~0(u6n3KtcZIKmcg(SA$v#P4$+crkJbNvy$(Y*YjOLZmF5SSkopNmQe*g~Xfc+KE zCk1+T2w(2pUU8b+=mX;Kb%@X|@1>Br}v3Kr9_<%bp1n8qE*x zh>P*vz_wWT_j z>6@*mO+KP3hUba8rL;L+v4M2JuU!dUDcz33=zpE|+`~+bIx-gj?lf< z-ua|IqEpJq(h}|3wYpGFujddy-rQNi9agWW} zHRwO-msPB3V_9lPO&Yo$T^m68GXOKJQOrr$k4#r}@FK8L-yj2rK)yE8~0_$jP|f7L-Z* zh)N%zt?3)yIp9}BxRJoxmRSPR*)NnmHEjG*j?69_n9}ll2t*+#oNflh3urNUNc~N9 zajG*j?eZ@r(qM)GG#<^uy@vJ$4sZzWVj{FSYWO6d%TKmFd@<2Z|7VHGSVX}B!1g%l zK&Z6+88F_IgQ}k|^p;s?T6RK*>cC_H5_^{u3=e{tFg1dpNQkdO~+o9}u!knOtHA$$4Kcu*>%%yQs7v<9p-c9I9=bol?% ct=igNQcg_rD6WQP6$1nMqOneu*45De0Y}m-3jhEB literal 0 HcmV?d00001 diff --git a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp index 41f6d2ca1..0f7f659b3 100644 --- a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp +++ b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp @@ -55,7 +55,31 @@ namespace GridKit return success.report(__func__); } - TestOutcome deadband() + TestOutcome deadband1() + { + TestStatus success = true; + + const ScalarT lower = scalar(-0.05); + const ScalarT upper = scalar(0.10); + const ScalarT tol = scalar(kSmoothTolerance); + + // Type 1 (no-offset) deadband returns ~0 inside the band ... + const ScalarT inside_input = scalar(0.02); + success *= (std::abs(Math::deadband1(inside_input, lower, upper)) < tol * tol); + + // ... and passes the input through unchanged (no offset) outside the band. + const ScalarT far_above_input = scalar(4.0); + const ScalarT far_below_input = scalar(-4.0); + + success *= std::isfinite(Math::deadband1(far_above_input, lower, upper)); + success *= within(Math::deadband1(far_above_input, lower, upper), far_above_input, tol); + success *= std::isfinite(Math::deadband1(far_below_input, lower, upper)); + success *= within(Math::deadband1(far_below_input, lower, upper), far_below_input, tol); + + return success.report(__func__); + } + + TestOutcome deadband2() { TestStatus success = true; @@ -69,12 +93,12 @@ namespace GridKit const ScalarT expected_below = below_input - lower; const ScalarT expected_above = above_input - upper; - success *= within(Math::deadband(below_input, lower, upper), expected_below, tol); - success *= (std::abs(Math::deadband(inside_input, lower, upper)) < tol * tol); - success *= within(Math::deadband(above_input, lower, upper), expected_above, tol); + success *= within(Math::deadband2(below_input, lower, upper), expected_below, tol); + success *= (std::abs(Math::deadband2(inside_input, lower, upper)) < tol * tol); + success *= within(Math::deadband2(above_input, lower, upper), expected_above, tol); - const ScalarT lower_breakpoint = Math::deadband(lower, lower, upper); - const ScalarT upper_breakpoint = Math::deadband(upper, lower, upper); + const ScalarT lower_breakpoint = Math::deadband2(lower, lower, upper); + const ScalarT upper_breakpoint = Math::deadband2(upper, lower, upper); success *= (lower_breakpoint < scalar(0.0)); success *= (upper_breakpoint > scalar(0.0)); @@ -82,7 +106,7 @@ namespace GridKit success *= (std::abs(upper_breakpoint) < tol); const ScalarT x = scalar(-0.4); - success *= (std::abs(Math::deadband(x, lower, upper) + success *= (std::abs(Math::deadband2(x, lower, upper) - (x - Math::clamp(x, lower, upper))) < scalar(kRoundoffTolerance)); @@ -91,17 +115,17 @@ namespace GridKit const ScalarT expected_far_above = far_above_input - upper; const ScalarT expected_far_below = far_below_input - lower; - success *= std::isfinite(Math::deadband(far_above_input, lower, upper)); - success *= within(Math::deadband(far_above_input, lower, upper), expected_far_above, tol); - success *= std::isfinite(Math::deadband(far_below_input, lower, upper)); - success *= within(Math::deadband(far_below_input, lower, upper), expected_far_below, tol); + success *= std::isfinite(Math::deadband2(far_above_input, lower, upper)); + success *= within(Math::deadband2(far_above_input, lower, upper), expected_far_above, tol); + success *= std::isfinite(Math::deadband2(far_below_input, lower, upper)); + success *= within(Math::deadband2(far_below_input, lower, upper), expected_far_below, tol); const ScalarT point = scalar(0.25); const ScalarT above_point = scalar(0.75); const ScalarT below_point = scalar(-0.25); - success *= (std::abs(Math::deadband(above_point, point, point) - (above_point - point)) + success *= (std::abs(Math::deadband2(above_point, point, point) - (above_point - point)) < scalar(kRoundoffTolerance)); - success *= (std::abs(Math::deadband(below_point, point, point) - (below_point - point)) + success *= (std::abs(Math::deadband2(below_point, point, point) - (below_point - point)) < scalar(kRoundoffTolerance)); return success.report(__func__); diff --git a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp index 233c36ddd..aecdc7ae8 100644 --- a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp +++ b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp @@ -7,7 +7,8 @@ int main() GridKit::Testing::SmoothnessIndicatorTests test; result += test.clamp(); - result += test.deadband(); + result += test.deadband1(); + result += test.deadband2(); result += test.limitIndicators(); result += test.slew(); result += test.linseg(); diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 932cb13b8..de1d7f639 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -78,6 +78,14 @@ target_link_libraries( GridKit::phasor_dynamics_bus_dependency_tracking GridKit::testing) +add_executable(test_phasor_governor_hygov runGovernorHygovTests.cpp) +target_link_libraries( + test_phasor_governor_hygov + GridKit::definitions + GridKit::phasor_dynamics_components + GridKit::phasor_dynamics_components_dependency_tracking + GridKit::testing) + add_executable(test_phasor_exciter_ieeet1 runExciterIeeet1Tests.cpp) target_link_libraries( test_phasor_exciter_ieeet1 @@ -135,6 +143,7 @@ add_test(NAME PhasorDynamicsBusToSignalAdapterTest COMMAND test_phasor_bustosign add_test(NAME PhasorDynamicsBranchTest COMMAND test_phasor_branch) add_test(NAME PhasorDynamicsGenrouTest COMMAND test_phasor_genrou) add_test(NAME PhasorDynamicsGovernorTgov1Test COMMAND test_phasor_governor_tgov1) +add_test(NAME PhasorDynamicsGovernorHygovTest COMMAND test_phasor_governor_hygov) 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) @@ -156,6 +165,7 @@ install( test_phasor_loadzip test_phasor_genrou test_phasor_governor_tgov1 + test_phasor_governor_hygov test_phasor_exciter_ieeet1 test_phasor_gensal test_phasor_exciter_sexspti diff --git a/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp b/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp new file mode 100644 index 000000000..d5495e229 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp @@ -0,0 +1,363 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class GovernorHygovTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + + static constexpr ScalarT kTol = static_cast(1.0e-8); + + TestOutcome constructionAndValidation() + { + TestStatus success = true; + + PhasorDynamics::Governor::Hygov hygov(makeHygovData()); + success *= (hygov.size() + == static_cast( + PhasorDynamics::Governor::HygovInternalVariables::MAXIMUM)); + success *= (hygov.getMonitor() != nullptr); + success *= (hygov.verify() == 0); + + auto source_default_hygov = makeHygovSourceDefaultData(); + PhasorDynamics::Governor::Hygov source_default_hygov_model( + source_default_hygov); + success *= (source_default_hygov_model.verify() == 0); + + auto bad_hygov = makeHygovData(); + bad_hygov.parameters[PhasorDynamics::Governor::HygovParameters::db2] = + static_cast(0.01); + PhasorDynamics::Governor::Hygov bad_hygov_model(bad_hygov); + success *= (bad_hygov_model.verify() > 0); + + auto bad_hygov_leadlag = makeHygovData(); + bad_hygov_leadlag.parameters[PhasorDynamics::Governor::HygovParameters::Tn] = + static_cast(0.1); + bad_hygov_leadlag.parameters[PhasorDynamics::Governor::HygovParameters::Tnp] = + static_cast(0.0); + PhasorDynamics::Governor::Hygov bad_hygov_leadlag_model( + bad_hygov_leadlag); + success *= (bad_hygov_leadlag_model.verify() > 0); + + auto bad_hygov_curve = makeHygovData(); + bad_hygov_curve.parameters[PhasorDynamics::Governor::HygovParameters::Gv2] = + static_cast(0.2); + PhasorDynamics::Governor::Hygov bad_hygov_curve_model(bad_hygov_curve); + success *= (bad_hygov_curve_model.verify() > 0); + + auto bad_hygov_trate = makeHygovData(); + bad_hygov_trate.parameters[PhasorDynamics::Governor::HygovParameters::Trate] = + static_cast(0.0); + PhasorDynamics::Governor::Hygov bad_hygov_trate_model(bad_hygov_trate); + success *= (bad_hygov_trate_model.verify() > 0); + + auto bad_hygov_base = makeHygovData(); + bad_hygov_base.parameters[PhasorDynamics::Governor::HygovParameters::mva] = + static_cast(0.0); + PhasorDynamics::Governor::Hygov bad_hygov_base_model(bad_hygov_base); + success *= (bad_hygov_base_model.verify() > 0); + + return success.report(__func__); + } + + TestOutcome signals() + { + using Var = PhasorDynamics::Governor::HygovInternalVariables; + + TestStatus success = true; + + PhasorDynamics::SignalNode pmech_node; + PhasorDynamics::SignalNode omega_node; + ScalarT pmech_value{0.40}; + ScalarT omega_value{0.0}; + IdxT pmech_index = 4; + IdxT omega_index = 5; + pmech_node.set(&pmech_value, &pmech_index); + omega_node.set(&omega_value, &omega_index); + + PhasorDynamics::Governor::Hygov hygov(makeHygovData()); + auto& hygov_signals = hygov.getSignals(); + hygov_signals + .template assignSignalNode( + &pmech_node); + hygov_signals + .template attachSignalNode( + &omega_node); + + success *= (hygov.allocate() == 0); + hygov.y()[index(Var::PMECH)] = pmech_value; + success *= (hygov.verify() == 0); + success *= (hygov.initialize() == 0); + success *= (hygov.tagDifferentiable() == 0); + success *= (hygov.evaluateResidual() == 0); + + success *= isEqual(hygov.y()[index(Var::Q)], static_cast(0.50), kTol); + success *= isEqual(hygov.y()[index(Var::G)], static_cast(0.50), kTol); + success *= isEqual(hygov.y()[index(Var::C)], static_cast(0.50), kTol); + success *= (hygov.tag()[index(Var::G)] == false); + + for (size_t i = 0; i < hygov.getResidual().size(); ++i) + { + success *= isEqual(hygov.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(hygov.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome sourceDefault() + { + using Var = PhasorDynamics::Governor::HygovInternalVariables; + + TestStatus success = true; + + PhasorDynamics::SignalNode pmech_node; + ScalarT pmech_value{0.40}; + IdxT pmech_index = 4; + pmech_node.set(&pmech_value, &pmech_index); + + PhasorDynamics::Governor::Hygov hygov(makeHygovSourceDefaultData()); + hygov.getSignals() + .template assignSignalNode( + &pmech_node); + + success *= (hygov.allocate() == 0); + hygov.y()[index(Var::PMECH)] = pmech_value; + success *= (hygov.verify() == 0); + success *= (hygov.initialize() == 0); + success *= (hygov.tagDifferentiable() == 0); + success *= (hygov.evaluateResidual() == 0); + + success *= isEqual(hygov.y()[index(Var::Q)], static_cast(0.50), kTol); + success *= isEqual(hygov.y()[index(Var::PGV)], static_cast(0.50), kTol); + success *= isEqual(hygov.y()[index(Var::G)], static_cast(0.50), kTol); + success *= (hygov.tag()[index(Var::XN)] == false); + + for (size_t i = 0; i < hygov.getResidual().size(); ++i) + { + success *= isEqual(hygov.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(hygov.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome baseConversion() + { + using Var = PhasorDynamics::Governor::HygovInternalVariables; + using Params = PhasorDynamics::Governor::HygovParameters; + + TestStatus success = true; + + PhasorDynamics::SignalNode pmech_node; + ScalarT pmech_value{0.40}; + IdxT pmech_index = 6; + pmech_node.set(&pmech_value, &pmech_index); + + auto data = makeHygovData(); + data.parameters[Params::Trate] = static_cast(50.0); + + PhasorDynamics::Governor::Hygov hygov(data); + hygov.getSignals().template assignSignalNode(&pmech_node); + + success *= (hygov.allocate() == 0); + hygov.y()[index(Var::PMECH)] = pmech_value; + success *= (hygov.verify() == 0); + success *= (hygov.initialize() == 0); + success *= (hygov.evaluateResidual() == 0); + + success *= isEqual(hygov.y()[index(Var::Q)], static_cast(0.90), kTol); + success *= isEqual(hygov.y()[index(Var::G)], static_cast(0.90), kTol); + success *= isEqual(hygov.y()[index(Var::PMECH)], static_cast(0.40), kTol); + success *= isEqual(pmech_node.read(), static_cast(0.40), kTol); + + for (size_t i = 0; i < hygov.getResidual().size(); ++i) + { + success *= isEqual(hygov.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(hygov.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome jsonParseAndSystemAssembly() + { + TestStatus success = true; + + std::istringstream input(R"json( +{ + "header": { + "format_version": 0, + "format_revision": 1, + "case_name": "hydro governor", + "case_description": "HYGOV 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": 10, "name": "Pmech" } + ], + "devices": [ + { + "class": "Genrou", + "ports": { "bus": 1, "pmech": 10 }, + "id": "GEN1", + "params": { + "p0": 0.3, "q0": 0.0, "H": 3.0, "D": 0.0, "Ra": 0.0, + "Tdop": 7.0, "Tdopp": 0.04, "Tqop": 0.75, "Tqopp": 0.05, + "Xd": 2.1, "Xdp": 0.2, "Xdpp": 0.18, "Xq": 0.5, "Xqp": 0.5, + "Xqpp": 0.18, "Xl": 0.15, "S10": 0.0, "S12": 0.0, + "mva": 100.0 + } + }, + { + "class": "Hygov", + "ports": { "pmech": 10 }, + "id": "HYG1", + "params": { + "Trate": 50.0, "mva": 100.0, "Rperm": 0.05, "Rtemp": 0.4, "Tr": 5.0, + "Tf": 0.2, "Tg": 0.0, "Velm": 0.5, "Gmax": 1.0, "Gmin": 0.0, + "Tw": 1.0, "At": 1.0, "Dturb": 0.0, "Qnl": 0.1, + "Tn": 0.0, "Tnp": 1.0, "db1": 0.0, "db2": 0.0, "Hdam": 1.0, + "Gv0": 0.0, "Gv1": 0.2, "Gv2": 0.4, "Gv3": 0.6, "Gv4": 0.8, "Gv5": 1.0, + "Pgv0": 0.0, "Pgv1": 0.2, "Pgv2": 0.4, "Pgv3": 0.6, "Pgv4": 0.8, "Pgv5": 1.0 + } + } + ] +} +)json"); + + auto data = PhasorDynamics::parseSystemModelData(input); + success *= (data.hygov.size() == 1); + const auto trate_param = + data.hygov[0].parameters.at(PhasorDynamics::Governor::HygovParameters::Trate); + const auto mva_param = + data.hygov[0].parameters.at(PhasorDynamics::Governor::HygovParameters::mva); + success *= (std::get(trate_param) == static_cast(50.0)); + success *= (std::get(mva_param) == static_cast(100.0)); + PhasorDynamics::SystemModel system(data); + success *= (system.allocate() == 0); + success *= (system.initialize() == 0); + const auto hygov_size = + static_cast(PhasorDynamics::Governor::HygovInternalVariables::MAXIMUM); + const auto hygov_offset = static_cast(system.size() - hygov_size); + success *= isEqual( + system.y()[hygov_offset + index(PhasorDynamics::Governor::HygovInternalVariables::Q)], + static_cast(0.70), + kTol); + success *= isEqual( + system.y()[hygov_offset + index(PhasorDynamics::Governor::HygovInternalVariables::G)], + static_cast(0.70), + kTol); + success *= (system.evaluateResidual() == 0); + success *= (system.size() == 35); + + return success.report(__func__); + } + + private: + static size_t index(PhasorDynamics::Governor::HygovInternalVariables variable) + { + return static_cast(variable); + } + + auto makeHygovData() -> PhasorDynamics::Governor::HygovData + { + using Params = PhasorDynamics::Governor::HygovParameters; + using Mon = PhasorDynamics::Governor::HygovMonitorableVariables; + + PhasorDynamics::Governor::HygovData data; + data.device_class = "Hygov"; + data.disambiguation_string = "hygov_test"; + data.monitored_variables.insert(Mon::pmech); + data.monitored_variables.insert(Mon::gate); + + data.parameters[Params::Trate] = static_cast(100.0); + data.parameters[Params::mva] = static_cast(100.0); + data.parameters[Params::Rperm] = static_cast(0.05); + data.parameters[Params::Rtemp] = static_cast(0.4); + data.parameters[Params::Tr] = static_cast(5.0); + data.parameters[Params::Tf] = static_cast(0.2); + data.parameters[Params::Tg] = static_cast(0.0); + data.parameters[Params::Velm] = static_cast(0.5); + data.parameters[Params::Gmax] = static_cast(1.0); + data.parameters[Params::Gmin] = static_cast(0.0); + data.parameters[Params::Tw] = static_cast(1.0); + data.parameters[Params::At] = static_cast(1.0); + data.parameters[Params::Dturb] = static_cast(0.0); + data.parameters[Params::Qnl] = static_cast(0.1); + data.parameters[Params::Tn] = static_cast(0.0); + data.parameters[Params::Tnp] = static_cast(1.0); + data.parameters[Params::db1] = static_cast(0.0); + data.parameters[Params::db2] = static_cast(0.0); + data.parameters[Params::Hdam] = static_cast(1.0); + data.parameters[Params::Gv0] = static_cast(0.0); + data.parameters[Params::Gv1] = static_cast(0.2); + data.parameters[Params::Gv2] = static_cast(0.4); + data.parameters[Params::Gv3] = static_cast(0.6); + data.parameters[Params::Gv4] = static_cast(0.8); + data.parameters[Params::Gv5] = static_cast(1.0); + data.parameters[Params::Pgv0] = static_cast(0.0); + data.parameters[Params::Pgv1] = static_cast(0.2); + data.parameters[Params::Pgv2] = static_cast(0.4); + data.parameters[Params::Pgv3] = static_cast(0.6); + data.parameters[Params::Pgv4] = static_cast(0.8); + data.parameters[Params::Pgv5] = static_cast(1.0); + + return data; + } + + auto makeHygovSourceDefaultData() -> PhasorDynamics::Governor::HygovData + { + using Params = PhasorDynamics::Governor::HygovParameters; + + auto data = makeHygovData(); + data.parameters[Params::Tn] = static_cast(0.0); + data.parameters[Params::Tnp] = static_cast(0.0); + data.parameters[Params::Gv0] = static_cast(0.0); + data.parameters[Params::Gv1] = static_cast(0.0); + data.parameters[Params::Gv2] = static_cast(0.0); + data.parameters[Params::Gv3] = static_cast(0.0); + data.parameters[Params::Gv4] = static_cast(0.0); + data.parameters[Params::Gv5] = static_cast(0.0); + data.parameters[Params::Pgv0] = static_cast(0.0); + data.parameters[Params::Pgv1] = static_cast(0.0); + data.parameters[Params::Pgv2] = static_cast(0.0); + data.parameters[Params::Pgv3] = static_cast(0.0); + data.parameters[Params::Pgv4] = static_cast(0.0); + data.parameters[Params::Pgv5] = static_cast(0.0); + + return data; + } + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp b/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp new file mode 100644 index 000000000..45154c841 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp @@ -0,0 +1,16 @@ +#include "GovernorHygovTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::GovernorHygovTests test; + + result += test.constructionAndValidation(); + result += test.signals(); + result += test.sourceDefault(); + result += test.baseConversion(); + result += test.jsonParseAndSystemAssembly(); + + return result.summary(); +}