diff --git a/docs/developer-guide/andes-compatibility.md b/docs/developer-guide/andes-compatibility.md
index 191785467..6f069cbc3 100644
--- a/docs/developer-guide/andes-compatibility.md
+++ b/docs/developer-guide/andes-compatibility.md
@@ -84,7 +84,8 @@ adapters still have attachment, parameter, equation, and validation gaps.
 | `IEEEG1`                                                                  | `IEEEG1`                | `GovernorIeeeG1`                                                              | **Implemented.** Frozen-ANDES equations and DYR order, initialization and perturbed equations, valve rate/position limits, analytic Jacobian, one- and two-machine attachment, and mixed synchronous-machine coupling are covered. Native ANDES import, alphanumeric DYR IDs, unequal machine bases, and a captured disturbed trajectory remain open.                                          |
 | `GAST`                                                                    | `GAST`                  | None exact                                                                    | **No direct analogue.** Add a gas-turbine governor implementation and DYR adapter.                                                                                                                                                                                                                                                                                                             |
 | `GGOV1`                                                                   | `TGOV1`                 | `GovernorTgov1`                                                               | **Planned compatibility approximation.** ANDES currently retains only `R` when converting this record. If reproduced, emit a diagnostic and keep exact GGOV1 support as a separate task.                                                                                                                                                                                                       |
-| `IEEEST`, `ST2CUT`                                                        | Same named ANDES models | Base `Stabilizer` only                                                        | **No direct analogue.** The base class is not a usable PSS implementation; add the two models and their generator/exciter signal connections.                                                                                                                                                                                                                                                  |
+| `ST2CUT`                                                                  | `ST2CUT`                | `StabilizerST2CUT`                                                            | **Partial.** Exact frozen-ANDES local modes 0/1/3/4/5, dual transducers, washout/lag, three lead-lag stages, output limits, voltage gating, DYR mapping, generator/exciter coupling, and dynamic load-step coverage are implemented. Remote `BUSR` inputs and ANDES modes 2/6 require cross-bus/frequency-derivative measurement routing and are rejected rather than approximated.            |
+| `IEEEST`                                                                  | `IEEEST`                | Base `Stabilizer` only                                                        | **No direct analogue.** The base class is not a usable PSS implementation; add the model and its generator/exciter signal connections.                                                                                                                                                                                                                                                         |
 | `REGCA1`, `REECA1`, `REECB1`, `REPCA1`                                    | Same named ANDES models | `GenModelInverter` is not equivalent                                          | **No direct analogue.** Add the coordinated renewable generator, electrical-control, and plant-control chain rather than flattening these records into the generic inverter.                                                                                                                                                                                                                   |
 | `WTDTA1`, `WTARA1`, `WTPTA1`, `WTTQA1`                                    | Same named ANDES models | None                                                                          | **No direct analogue.** Add drive-train, aerodynamic, pitch, and torque-control submodels with their shared interfaces.                                                                                                                                                                                                                                                                        |
 | `Toggle`, `Fault`                                                         | ANDES event models      | GridDyn event/action and fault mechanisms                                     | **Partial conceptually.** Define DYR record schemas and translate target resolution, timing, status changes, and fault clearing semantics; no adapter exists.                                                                                                                                                                                                                                  |
@@ -170,14 +171,14 @@ seconds, and continuing through 2.0 seconds. The five controller chains are:
 
 ### PR status and dependency order
 
-| PR  | Deliverable                                                                                                   | Status                                         | Depends on |
-| --- | ------------------------------------------------------------------------------------------------------------- | ---------------------------------------------- | ---------- |
-| 1   | Complete GENROU equations, saturation utility, DYR mapping, and five-machine initialization reference         | Current changes; treat as complete after merge | none       |
-| 2   | Robust RAW/DYR identity and base handling, IEEE 14 power-flow parity, and a controller-free GENROU trajectory | Planned                                        | PR 1       |
-| 3   | Generator/controller signal plumbing, PSS-to-exciter routing, and validated `TGOV1`                           | Implemented                                    | PR 2       |
-| 4   | Complete `ESST3A` and `EXST1` models plus DYR adapters                                                        | Models/adapters implemented; trajectories open | PR 3       |
-| 5   | Complete `IEEEG1`, `ST2CUT`, and `IEEEST` models plus DYR adapters                                            | `IEEEG1` implemented; stabilizers planned      | PR 4       |
-| 6   | DYR `Toggle`, complete IEEE 14 initialization/equilibrium, and the two-second trajectory regression           | Planned                                        | PR 5       |
+| PR  | Deliverable                                                                                                   | Status                                                   | Depends on |
+| --- | ------------------------------------------------------------------------------------------------------------- | -------------------------------------------------------- | ---------- |
+| 1   | Complete GENROU equations, saturation utility, DYR mapping, and five-machine initialization reference         | Current changes; treat as complete after merge           | none       |
+| 2   | Robust RAW/DYR identity and base handling, IEEE 14 power-flow parity, and a controller-free GENROU trajectory | Planned                                                  | PR 1       |
+| 3   | Generator/controller signal plumbing, PSS-to-exciter routing, and validated `TGOV1`                           | Implemented                                              | PR 2       |
+| 4   | Complete `ESST3A` and `EXST1` models plus DYR adapters                                                        | Models/adapters implemented; trajectories open           | PR 3       |
+| 5   | Complete `IEEEG1`, `ST2CUT`, and `IEEEST` models plus DYR adapters                                            | `IEEEG1` implemented; `ST2CUT` partial; `IEEEST` planned | PR 4       |
+| 6   | DYR `Toggle`, complete IEEE 14 initialization/equilibrium, and the two-second trajectory regression           | Planned                                                  | PR 5       |
 
 PRs must merge in this order. Each PR must pass without relying on production
 code from a later PR. Minimized fixtures and test-only probe models are
@@ -479,6 +480,33 @@ equations, Jacobian, limit/root transitions, and a minimized coupled ANDES
 trajectory. Add explicit remote-bus and machine/system-base tests where the
 schema permits them.
 
+#### ST2CUT implementation and parameter contract
+
+`StabilizerST2CUT` follows frozen ANDES v2.0.0
+`andes/models/pss/st2cut.py`; its DYR adapter uses the exact frozen
+`psse-dyr.yaml` order `BUS, ID, MODE, BUSR, MODE2, BUSR2, K1, K2, T1, T2,
+T3, T4, T5, T6, T7, T8, T9, T10, LSMAX, LSMIN, VCU, VCL`. The two input
+transducers feed the ANDES washout-or-lag block, three lead-lag blocks, the
+`[LSMIN, LSMAX]` output limiter, and the initialized-voltage-relative
+`[VCL, VCU]` output gate. As in ANDES, zero `VCU`/`VCL` is mapped to
+`+999`/`-999` for an effectively disabled voltage gate.
+
+GridDyn currently provides local rotor speed, terminal voltage, mechanical
+power, and electrical power. Therefore local modes 0, 1, 3, 4, and 5 are
+implemented and tested. ANDES mode 2 (remote bus-frequency), mode 6 (terminal
+voltage derivative), and nonzero `BUSR`/`BUSR2` need a generic cross-bus
+measurement-routing interface that GridDyn does not yet have; these records
+are rejected at parsing/initialization rather than silently using a local
+signal. `T3=0` uses the exact ANDES lag form. The nonzero lag denominators
+`T1`, `T2`, `T4`, `T6`, `T8`, and `T10` are required; zero-denominator block
+bypasses remain unsupported.
+
+Focused tests cover DYR mapping, initialization and perturbed equations,
+analytic-Jacobian finite differences, limiter and voltage-gate root
+transitions, GENROU plus exciter coupling, and a load-step dynamic run whose
+ST2CUT states must evolve. A captured ANDES trajectory and support for the
+remote/derivative measurement modes remain open.
+
 Merge gate: a no-event IEEE 14 controller subset loads all five governors and
 all three stabilizers with no unsupported-record diagnostics; each controller
 also passes its independent minimized numerical comparison. Do not add the
diff --git a/src/fileInput/gridDynReadDYR.cpp b/src/fileInput/gridDynReadDYR.cpp
index 01283fcea..52a516133 100644
--- a/src/fileInput/gridDynReadDYR.cpp
+++ b/src/fileInput/gridDynReadDYR.cpp
@@ -17,8 +17,10 @@
 #include "griddyn/Generator.h"
 #include "griddyn/Governor.h"
 #include "griddyn/GridBus.h"
+#include "griddyn/Stabilizer.h"
 #include "griddyn/generators/DynamicGenerator.h"
 #include "griddyn/governors/GovernorIeeeG1.h"
+#include "griddyn/stabilizers/StabilizerST2CUT.h"
 #include <cmath>
 #include <fstream>
 #include <iostream>
@@ -33,6 +35,7 @@ namespace {
     void loadEXST1(CoreObject* parentObject, stringVec& tokens);
     void loadTGOV1(CoreObject* parentObject, stringVec& tokens);
     void loadIEEEG1(CoreObject* parentObject, stringVec& tokens);
+    void loadST2CUT(CoreObject* parentObject, stringVec& tokens);
     void loadEXDC2(CoreObject* parentObject, stringVec& tokens);
     void loadSEXS(CoreObject* parentObject, stringVec& tokens);
 }  // namespace
@@ -87,6 +90,8 @@ void loadDyr(CoreObject* parentObject,
             loadTGOV1(parentObject, lineTokens);
         } else if (type == "'IEEEG1'") {
             loadIEEEG1(parentObject, lineTokens);
+        } else if (type == "'ST2CUT'") {
+            loadST2CUT(parentObject, lineTokens);
         } else if (type == "'SEXS'") {
             loadSEXS(parentObject, lineTokens);
         } else {
@@ -378,6 +383,47 @@ namespace {
                                                 governors::GovernorIeeeG1::lpOutput);
         }
     }
+
+    void loadST2CUT(CoreObject* parentObject, stringVec& tokens)
+    {
+        if (tokens.size() != 23U) {
+            throw InvalidParameterValue("ST2CUT DYR record must contain 23 fields");
+        }
+        const int busId = std::stoi(tokens[0]);
+        const auto* bus = static_cast<GridBus*>(parentObject->findByUserID("bus", busId));
+        const int genId = std::stoi(tokens[2]);
+        auto* generator = dynamic_cast<DynamicGenerator*>(bus->getGen(genId - 1));
+        if (generator == nullptr) {
+            throw InvalidParameterValue("ST2CUT requires a dynamic generator");
+        }
+
+        const auto params = gmlc::utilities::str2vector(tokens, kNullVal);
+        auto* stabilizer = new stabilizers::StabilizerST2CUT();
+        // Exact frozen ANDES psse-dyr.yaml order after BUS and ID:
+        // MODE, BUSR, MODE2, BUSR2, K1, K2, T1, T2, T3, T4,
+        // T5, T6, T7, T8, T9, T10, LSMAX, LSMIN, VCU, VCL.
+        stabilizer->set("mode", params[3]);
+        stabilizer->set("busr", params[4]);
+        stabilizer->set("mode2", params[5]);
+        stabilizer->set("busr2", params[6]);
+        stabilizer->set("k1", params[7]);
+        stabilizer->set("k2", params[8]);
+        stabilizer->set("t1", params[9]);
+        stabilizer->set("t2", params[10]);
+        stabilizer->set("t3", params[11]);
+        stabilizer->set("t4", params[12]);
+        stabilizer->set("t5", params[13]);
+        stabilizer->set("t6", params[14]);
+        stabilizer->set("t7", params[15]);
+        stabilizer->set("t8", params[16]);
+        stabilizer->set("t9", params[17]);
+        stabilizer->set("t10", params[18]);
+        stabilizer->set("lsmax", params[19]);
+        stabilizer->set("lsmin", params[20]);
+        stabilizer->set("vcu", params[21]);
+        stabilizer->set("vcl", params[22]);
+        generator->add(stabilizer);
+    }
 }  // namespace
 
 }  // namespace griddyn
diff --git a/src/griddyn/CMakeLists.txt b/src/griddyn/CMakeLists.txt
index b04b0b188..d2839003d 100644
--- a/src/griddyn/CMakeLists.txt
+++ b/src/griddyn/CMakeLists.txt
@@ -81,6 +81,7 @@ set(gridSource_sources
 set(exciter_headers
     Exciter.h
     Stabilizer.h
+    stabilizers/StabilizerST2CUT.h
     exciters/ExciterDC1A.h
     exciters/ExciterDC2A.h
     exciters/ExciterESST3A.h
@@ -101,6 +102,7 @@ set(exciter_sources
     exciters/ExciterIEEEtype2.cpp
     exciters/ExciterSEXS.cpp
     stabilizers/Stabilizer.cpp
+    stabilizers/StabilizerST2CUT.cpp
 )
 
 set(gov_headers
diff --git a/src/griddyn/Stabilizer.h b/src/griddyn/Stabilizer.h
index 1035fa06c..780a80b51 100644
--- a/src/griddyn/Stabilizer.h
+++ b/src/griddyn/Stabilizer.h
@@ -8,8 +8,15 @@
 
 #include "GridSubModel.h"
 #include <string>
+#include <vector>
 
 namespace griddyn {
+inline constexpr int pssOmegaInLocation = 0;
+inline constexpr int pssVoltageInLocation = 1;
+inline constexpr int pssPmechInLocation = 2;
+inline constexpr int pssElectricalPowerInLocation = 3;
+inline constexpr count_t pssInputCount = 4;
+
 class Stabilizer: public GridSubModel {
   public:
   protected:
@@ -48,6 +55,9 @@ class Stabilizer: public GridSubModel {
                             double deriv[],
                             const SolverMode& sMode) override;
 
+    virtual const std::vector<stringVec>& inputNames() const override;
+    virtual const std::vector<stringVec>& outputNames() const override;
+
     virtual index_t findIndex(std::string_view field, const SolverMode& sMode) const override;
 };
 
diff --git a/src/griddyn/generators/DynamicGenerator.cpp b/src/griddyn/generators/DynamicGenerator.cpp
index 4b9ab5c5c..3d2b933f5 100644
--- a/src/griddyn/generators/DynamicGenerator.cpp
+++ b/src/griddyn/generators/DynamicGenerator.cpp
@@ -97,6 +97,19 @@ namespace {
         }
     }
 
+    void copyPssInputs(const IOdata& machineSignals,
+                       double omega,
+                       double voltage,
+                       double mechanicalPower,
+                       IOdata& pssInputs)
+    {
+        pssInputs[pssOmegaInLocation] = omega;
+        pssInputs[pssVoltageInLocation] = voltage;
+        pssInputs[pssPmechInLocation] = mechanicalPower;
+        pssInputs[pssElectricalPowerInLocation] =
+            machineSignals[static_cast<index_t>(MachineControllerSignal::ELECTRICAL_TORQUE)];
+    }
+
     const auto& getDynModelFromStringMap()
     {
         static const std::map<std::string_view, DynamicGenerator::DynModel, std::less<>>
@@ -325,17 +338,23 @@ void DynamicGenerator::dynObjectInitializeB(const IOdata& inputs,
     }
 
     if ((pss != nullptr) && (pss->isEnabled())) {
-        modelInputs[0] = systemBaseFrequency;
-        modelInputs[1] = kNullVal;
+        IOdata pssInputs(pssInputCount, 0.0);
+        pssInputs[pssOmegaInLocation] = 1.0;
+        pssInputs[pssVoltageInLocation] = voltage;
+        pssInputs[pssPmechInLocation] = m_Pmech;
+        pssInputs[pssElectricalPowerInLocation] = m_Pmech;
         localDesiredOutput[0] = 0;
-        pss->dynInitializeB(modelInputs, desiredOutput, computedFieldSet);
+        pss->dynInitializeB(pssInputs, localDesiredOutput, computedFieldSet);
         //    pss->guessState (prevTime, m_state.data (), m_dstate_dt.data (), cLocalbSolverMode);
     }
 
     modelInputs.resize(0);
     localDesiredOutput.resize(0);
     for (auto* sub : getSubObjects()) {
-        if (sub->locIndex < 4) {
+        // The machine, exciter, governor, and PSS above require their own
+        // controller input contracts.  Do not initialize them again through
+        // the generic empty-input path.
+        if (sub->locIndex <= PSS_LOC) {
             continue;
         }
         if (sub->isEnabled()) {
@@ -623,7 +642,19 @@ void DynamicGenerator::timestep(CoreTime time, const IOdata& inputs, const Solve
         }
 
         if ((pss != nullptr) && (pss->isEnabled())) {
-            pss->timestep(time, inputs, sMode);
+            const IOdata genModelInputs{inputs[VOLTAGE_IN_LOCATION],
+                                        inputs[ANGLE_IN_LOCATION],
+                                        m_Eft,
+                                        m_Pmech};
+            IOdata pssInputs(pssInputCount, 0.0);
+            copyPssInputs(genModel->getMachineControllerSignals(genModelInputs,
+                                                                emptyStateData,
+                                                                cLocalSolverMode),
+                          omega,
+                          inputs[VOLTAGE_IN_LOCATION],
+                          m_Pmech,
+                          pssInputs);
+            pss->timestep(time, pssInputs, sMode);
         }
 
         if ((ext != nullptr) && (ext->isEnabled())) {
@@ -1023,7 +1054,7 @@ void DynamicGenerator::jacobianElements(const IOdata& inputs,
     // compute the Jacobian
     for (auto* sub : getSubObjects()) {
         if (sub->isEnabled()) {
-            if ((sub == ext) && (ext->numInputs() > exciterMachineSignalBase)) {
+            if (((sub == ext) && (ext->numInputs() > exciterMachineSignalBase)) || (sub == pss)) {
                 MatrixDataCustomWriteOnly<double> translatedMatrix;
                 translatedMatrix.setFunction(
                     [&matrixDataValue,
@@ -1182,6 +1213,7 @@ DynamicGenerator::SubModelInputs::SubModelInputs(): inputs(6)
     inputs[GEN_MODEL_LOC].resize(4);
     inputs[EXCITER_LOC].resize(exciterInputCount);
     inputs[GOVERNOR_LOC].resize(3);
+    inputs[PSS_LOC].resize(pssInputCount);
 }
 
 DynamicGenerator::SubModelInputLocs::SubModelInputLocs():
@@ -1190,6 +1222,7 @@ DynamicGenerator::SubModelInputLocs::SubModelInputLocs():
     inputLocs[GEN_MODEL_LOC].resize(4);
     inputLocs[EXCITER_LOC].resize(exciterInputCount);
     inputLocs[GOVERNOR_LOC].resize(3);
+    inputLocs[PSS_LOC].resize(pssInputCount);
 
     genModelInputLocsExternal[genModelEftInLocation] = kNullLocation;
     genModelInputLocsExternal[genModelPmechInLocation] = kNullLocation;
@@ -1254,10 +1287,18 @@ void DynamicGenerator::generateSubModelInputs(const IOdata& inputs,
     subInputs.inputs[EXCITER_LOC][exciterPmechInLocation] = pmech;
     subInputs.inputs[EXCITER_LOC][exciterVsetInLocation] =
         vSetControlUpdate(inputs, stateDataValue, sMode);
-    copyMachineSignals(genModel->getMachineControllerSignals(subInputs.inputs[GEN_MODEL_LOC],
-                                                             stateDataValue,
-                                                             sMode),
-                       subInputs.inputs[EXCITER_LOC]);
+    const auto machineSignals =
+        genModel->getMachineControllerSignals(subInputs.inputs[GEN_MODEL_LOC],
+                                              stateDataValue,
+                                              sMode);
+    copyMachineSignals(machineSignals, subInputs.inputs[EXCITER_LOC]);
+    if ((pss != nullptr) && (pss->isEnabled())) {
+        copyPssInputs(machineSignals,
+                      subInputs.inputs[GOVERNOR_LOC][govOmegaInLocation],
+                      subInputs.inputs[EXCITER_LOC][exciterVoltageInLocation],
+                      pmech,
+                      subInputs.inputs[PSS_LOC]);
+    }
     subInputs.inputs[EXCITER_LOC][exciterVssInLocation] = 0.0;
     if ((pss != nullptr) && (pss->isEnabled()) && (pss->numOutputs() > 0)) {
         subInputs.inputs[EXCITER_LOC][exciterVssInLocation] =
@@ -1330,6 +1371,17 @@ void DynamicGenerator::generateSubModelInputLocs(const IOlocs& inputLocs,
     subInputLocs.genModelInputLocsInternal[genModelPmechInLocation] =
         subInputLocs.inputLocs[GEN_MODEL_LOC][genModelPmechInLocation];
 
+    if ((pss != nullptr) && (pss->isEnabled())) {
+        index_t omegaLocation;
+        genModel->getFreq(stateDataValue, sMode, &omegaLocation);
+        subInputLocs.inputLocs[PSS_LOC][pssOmegaInLocation] = omegaLocation;
+        subInputLocs.inputLocs[PSS_LOC][pssVoltageInLocation] = inputLocs[VOLTAGE_IN_LOCATION];
+        subInputLocs.inputLocs[PSS_LOC][pssPmechInLocation] =
+            subInputLocs.inputLocs[GEN_MODEL_LOC][genModelPmechInLocation];
+        subInputLocs.inputLocs[PSS_LOC][pssElectricalPowerInLocation] =
+            machineSignalLocation(static_cast<index_t>(MachineControllerSignal::ELECTRICAL_TORQUE));
+    }
+
     if (isoc != nullptr) {
         subInputLocs.inputLocs[ISOC_CONTROL_LOC][0] =
             subInputLocs.inputLocs[GOVERNOR_LOC][govOmegaInLocation];
diff --git a/src/griddyn/stabilizers/Stabilizer.cpp b/src/griddyn/stabilizers/Stabilizer.cpp
index 5f13c63eb..367290193 100644
--- a/src/griddyn/stabilizers/Stabilizer.cpp
+++ b/src/griddyn/stabilizers/Stabilizer.cpp
@@ -8,18 +8,28 @@
 
 #include "../Generator.h"
 #include "../GridBus.h"
+#include "StabilizerST2CUT.h"
 #include "core/CoreObjectTemplates.hpp"
 #include "core/ObjectFactoryTemplates.hpp"
 #include <cmath>
 #include <string>
+#include <vector>
 
 namespace griddyn {
-static const TypeFactory<Stabilizer> gf("pss", std::to_array<std::string_view>({"basic"}));
+static const TypeFactory<Stabilizer> STABILIZER_FACTORY("pss",
+                                                        std::to_array<std::string_view>({"basic"}));
+static ChildTypeFactory<stabilizers::StabilizerST2CUT, Stabilizer>
+    gSt2cutFactory("pss", std::to_array<std::string_view>({"st2cut"}));
 
-Stabilizer::Stabilizer(const std::string& objName): GridSubModel(objName) {}
+Stabilizer::Stabilizer(const std::string& objName):
+    GridSubModel(objName), mp_Tw(0.0), mp_Teps(0.0), mp_Kw(0.0), mp_Kp(0.0), mp_Kv(0.0),
+    mp_Smax(0.0), mp_Smin(0.0)
+{
+    m_inputSize = pssInputCount;
+}
 CoreObject* Stabilizer::clone(CoreObject* obj) const
 {
-    auto pss = cloneBase<Stabilizer, GridSubModel>(this, obj);
+    auto* pss = cloneBase<Stabilizer, GridSubModel>(this, obj);
     if (pss == nullptr) {
         return obj;
     }
@@ -51,7 +61,7 @@ index_t Stabilizer::findIndex(std::string_view /*field*/, const SolverMode& /*sM
 
 void Stabilizer::set(std::string_view param, std::string_view val)
 {
-    return CoreObject::set(param, val);
+    CoreObject::set(param, val);
 }
 // set parameters
 void Stabilizer::set(std::string_view param, double val, units::unit unitType)
@@ -79,4 +89,19 @@ void Stabilizer::derivative(const IOdata& /*inputs*/,
 {
 }
 
+const std::vector<stringVec>& Stabilizer::inputNames() const
+{
+    static const std::vector<stringVec> inputNamesStr{{"omega", "frequency", "w"},
+                                                      {"voltage", "v", "volt"},
+                                                      {"pmech", "mechanicalpower"},
+                                                      {"pe", "electricalpower", "te"}};
+    return inputNamesStr;
+}
+
+const std::vector<stringVec>& Stabilizer::outputNames() const
+{
+    static const std::vector<stringVec> outputNamesStr{{"vss", "stabilizersignal"}};
+    return outputNamesStr;
+}
+
 }  // namespace griddyn
diff --git a/src/griddyn/stabilizers/StabilizerST2CUT.cpp b/src/griddyn/stabilizers/StabilizerST2CUT.cpp
new file mode 100644
index 000000000..95f7a42ff
--- /dev/null
+++ b/src/griddyn/stabilizers/StabilizerST2CUT.cpp
@@ -0,0 +1,622 @@
+/*
+ * Copyright (c) 2014-2026, Lawrence Livermore National Security
+ * See the top-level NOTICE for additional details. All rights reserved.
+ * SPDX-License-Identifier: BSD-3-Clause
+ */
+
+#include "StabilizerST2CUT.h"
+
+#include "core/CoreExceptions.h"
+#include "core/CoreObjectTemplates.hpp"
+#include "utilities/MatrixData.hpp"
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <string>
+#include <vector>
+
+namespace griddyn::stabilizers {
+namespace {
+    constexpr index_t transducer1State = 0;
+    constexpr index_t transducer2State = 1;
+    constexpr index_t washoutState = 2;
+    constexpr index_t leadLag1State = 3;
+    constexpr index_t leadLag2State = 4;
+    constexpr index_t leadLag3State = 5;
+}  // namespace
+
+StabilizerST2CUT::StabilizerST2CUT(const std::string& objName): Stabilizer(objName)
+{
+    m_inputSize = pssInputCount;
+    m_outputSize = 1;
+}
+
+CoreObject* StabilizerST2CUT::clone(CoreObject* obj) const
+{
+    auto* stabilizerClone = cloneBase<StabilizerST2CUT, Stabilizer>(this, obj);
+    if (stabilizerClone == nullptr) {
+        return obj;
+    }
+    stabilizerClone->mode1 = mode1;
+    stabilizerClone->remoteBus1 = remoteBus1;
+    stabilizerClone->mode2 = mode2;
+    stabilizerClone->remoteBus2 = remoteBus2;
+    stabilizerClone->K1 = K1;
+    stabilizerClone->K2 = K2;
+    stabilizerClone->T1 = T1;
+    stabilizerClone->T2 = T2;
+    stabilizerClone->T3 = T3;
+    stabilizerClone->T4 = T4;
+    stabilizerClone->T5 = T5;
+    stabilizerClone->T6 = T6;
+    stabilizerClone->T7 = T7;
+    stabilizerClone->T8 = T8;
+    stabilizerClone->T9 = T9;
+    stabilizerClone->T10 = T10;
+    stabilizerClone->Lsmax = Lsmax;
+    stabilizerClone->Lsmin = Lsmin;
+    stabilizerClone->Vcu = Vcu;
+    stabilizerClone->Vcl = Vcl;
+    return stabilizerClone;
+}
+
+bool StabilizerST2CUT::supportedMode(int mode)
+{
+    return (mode == 0) || (mode == 1) || (mode == 3) || (mode == 4) || (mode == 5);
+}
+
+void StabilizerST2CUT::dynObjectInitializeA(CoreTime /*time0*/, std::uint32_t /*flags*/)
+{
+    const std::array<double, 16> parameters{
+        K1, K2, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, Lsmax, Lsmin, Vcu, Vcl};
+    if (!std::all_of(parameters.begin(),
+                     parameters.end(),
+                     [](double value) { return std::isfinite(value); }) ||
+        !supportedMode(mode1) || !supportedMode(mode2) || (remoteBus1 != 0) || (remoteBus2 != 0) ||
+        (T1 <= 0.0) || (T2 <= 0.0) || (T3 < 0.0) || (T4 <= 0.0) || (T6 <= 0.0) || (T8 <= 0.0) ||
+        (T5 < 0.0) || (T7 < 0.0) || (T9 < 0.0) || (T10 <= 0.0) || (Lsmax < Lsmin) || (Vcu < Vcl)) {
+        throw InvalidParameterValue("ST2CUT modes, time constants, limits, or remote-bus inputs");
+    }
+    offsets.local().local.algSize = 1;
+    offsets.local().local.diffSize = 6;
+    offsets.local().local.algRoots = 4;
+    offsets.local().local.jacSize = 28;
+}
+
+void StabilizerST2CUT::dynObjectInitializeB(const IOdata& inputs,
+                                            const IOdata& /*desiredOutput*/,
+                                            IOdata& /*fieldSet*/)
+{
+    const auto validInput = [&inputs](index_t index) {
+        return std::isfinite(inputs[index]) && (std::abs(inputs[index]) < 1e20);
+    };
+    const auto selectedInputAvailable = [&validInput](int mode) {
+        switch (mode) {
+            case 0:
+                return true;
+            case 1:
+                return validInput(pssOmegaInLocation);
+            case 3:
+                return validInput(pssElectricalPowerInLocation);
+            case 4:
+                return validInput(pssPmechInLocation);
+            case 5:
+                return validInput(pssVoltageInLocation);
+            default:
+                return false;
+        }
+    };
+    if ((inputs.size() < pssInputCount) || !validInput(pssVoltageInLocation) ||
+        !selectedInputAvailable(mode1) || !selectedInputAvailable(mode2)) {
+        throw InvalidParameterValue("ST2CUT controller inputs");
+    }
+    initialVoltage = inputs[pssVoltageInLocation];
+    initialPmech = ((mode1 == 4) || (mode2 == 4)) ? inputs[pssPmechInLocation] : 0.0;
+    const double input1 = selectedInput(inputs, mode1);
+    const double input2 = selectedInput(inputs, mode2);
+    m_state[0] = 0.0;
+    double* state = m_state.data() + 1;
+    state[transducer1State] = K1 * input1;
+    state[transducer2State] = K2 * input2;
+    const double summedInput = state[transducer1State] + state[transducer2State];
+    state[washoutState] = summedInput;
+    state[leadLag1State] = 0.0;
+    state[leadLag2State] = 0.0;
+    state[leadLag3State] = 0.0;
+    std::fill(m_dstate_dt.begin(), m_dstate_dt.end(), 0.0);
+    updateLimitFlags(inputs, state);
+}
+
+double StabilizerST2CUT::selectedInput(const IOdata& inputs, int mode) const
+{
+    switch (mode) {
+        case 0:
+            return 0.0;
+        case 1:
+            return inputs[pssOmegaInLocation] - 1.0;
+        case 3:
+            return inputs[pssElectricalPowerInLocation];
+        case 4:
+            return inputs[pssPmechInLocation] - initialPmech;
+        case 5:
+            return inputs[pssVoltageInLocation];
+        default:
+            return 0.0;
+    }
+}
+
+double StabilizerST2CUT::washoutOutput(const double state[]) const
+{
+    const double summedInput = state[transducer1State] + state[transducer2State];
+    if (T3 <= 0.0) {
+        return state[washoutState];
+    }
+    return (T3 / T4) * (summedInput - state[washoutState]);
+}
+
+double StabilizerST2CUT::leadLagOutput(double input, double state, double leadTime, double lagTime)
+{
+    return state + ((leadTime / lagTime) * (input - state));
+}
+
+double StabilizerST2CUT::unlimitedOutput(const double state[]) const
+{
+    const double output1 = leadLagOutput(washoutOutput(state), state[leadLag1State], T5, T6);
+    const double output2 = leadLagOutput(output1, state[leadLag2State], T7, T8);
+    return leadLagOutput(output2, state[leadLag3State], T9, T10);
+}
+
+bool StabilizerST2CUT::voltageEnabled(const IOdata& inputs) const
+{
+    return (inputs[pssVoltageInLocation] >= initialVoltage + Vcl) &&
+        (inputs[pssVoltageInLocation] <= initialVoltage + Vcu);
+}
+
+int StabilizerST2CUT::outputLimitStatus(const double state[]) const
+{
+    const double outputValue = unlimitedOutput(state);
+    if (outputValue >= Lsmax) {
+        return 1;
+    }
+    return (outputValue <= Lsmin) ? -1 : 0;
+}
+
+double StabilizerST2CUT::output(const IOdata& inputs, const double state[]) const
+{
+    if (!voltageEnabled(inputs)) {
+        return 0.0;
+    }
+    return std::clamp(unlimitedOutput(state), Lsmin, Lsmax);
+}
+
+bool StabilizerST2CUT::updateLimitFlags(const IOdata& inputs, const double state[])
+{
+    const int outputStatus = outputLimitStatus(state);
+    const bool gated = !voltageEnabled(inputs);
+    const bool changed = (opFlags[OUTPUT_LIMITED] != (outputStatus != 0)) ||
+        (opFlags[OUTPUT_LIMIT_HIGH] != (outputStatus > 0)) || (opFlags[VOLTAGE_GATED] != gated);
+    opFlags.set(OUTPUT_LIMITED, outputStatus != 0);
+    opFlags.set(OUTPUT_LIMIT_HIGH, outputStatus > 0);
+    opFlags.set(VOLTAGE_GATED, gated);
+    return changed;
+}
+
+void StabilizerST2CUT::residual(const IOdata& inputs,
+                                const StateData& stateData,
+                                double resid[],
+                                const SolverMode& sMode)
+{
+    const auto locations = offsets.getLocations(stateData, resid, sMode, this);
+    if (hasAlgebraic(sMode)) {
+        locations.destLoc[0] = output(inputs, locations.diffStateLoc) - locations.algStateLoc[0];
+    }
+    if (!hasDifferential(sMode)) {
+        return;
+    }
+    derivative(inputs, stateData, resid, sMode);
+    for (index_t index = 0; index < locations.diffSize; ++index) {
+        locations.destDiffLoc[index] -= locations.dstateLoc[index];
+    }
+}
+
+void StabilizerST2CUT::derivative(const IOdata& inputs,
+                                  const StateData& stateData,
+                                  double deriv[],
+                                  const SolverMode& sMode)
+{
+    if (!hasDifferential(sMode)) {
+        return;
+    }
+    const auto locations = offsets.getLocations(stateData, deriv, sMode, this);
+    const double* state = locations.diffStateLoc;
+    double* stateDerivative = locations.destDiffLoc;
+    const double input1 = selectedInput(inputs, mode1);
+    const double input2 = selectedInput(inputs, mode2);
+    const double summedInput = state[transducer1State] + state[transducer2State];
+    const double washout = washoutOutput(state);
+    const double output1 = leadLagOutput(washout, state[leadLag1State], T5, T6);
+    const double output2 = leadLagOutput(output1, state[leadLag2State], T7, T8);
+
+    stateDerivative[transducer1State] = ((K1 * input1) - state[transducer1State]) / T1;
+    stateDerivative[transducer2State] = ((K2 * input2) - state[transducer2State]) / T2;
+    stateDerivative[washoutState] = (summedInput - state[washoutState]) / T4;
+    stateDerivative[leadLag1State] = (washout - state[leadLag1State]) / T6;
+    stateDerivative[leadLag2State] = (output1 - state[leadLag2State]) / T8;
+    stateDerivative[leadLag3State] = (output2 - state[leadLag3State]) / T10;
+}
+
+void StabilizerST2CUT::jacobianElements(const IOdata& inputs,
+                                        const StateData& stateData,
+                                        MatrixData<double>& matrixData,
+                                        const IOlocs& inputLocs,
+                                        const SolverMode& sMode)
+{
+    const auto locations = offsets.getLocations(stateData, sMode, this);
+    const index_t refAlg = locations.algOffset;
+    const index_t refDiff = locations.diffOffset;
+    const double* state = locations.diffStateLoc;
+
+    const double washoutScale = (T3 <= 0.0) ? 0.0 : T3 / T4;
+    const double output1Gain = T5 / T6;
+    const double output2Gain = T7 / T8;
+    const double output3Gain = T9 / T10;
+    const double ll1StateGain = 1.0 - output1Gain;
+    const double ll2StateGain = 1.0 - output2Gain;
+    const double ll3StateGain = 1.0 - output3Gain;
+
+    if (hasAlgebraic(sMode)) {
+        matrixData.assign(refAlg, refAlg, -1.0);
+        if (voltageEnabled(inputs) && (outputLimitStatus(state) == 0)) {
+            matrixData.assign(refAlg, refDiff + leadLag3State, ll3StateGain);
+            matrixData.assign(refAlg, refDiff + leadLag2State, output3Gain * ll2StateGain);
+            matrixData.assign(refAlg,
+                              refDiff + leadLag1State,
+                              output3Gain * output2Gain * ll1StateGain);
+            if (T3 <= 0.0) {
+                matrixData.assign(refAlg,
+                                  refDiff + washoutState,
+                                  output3Gain * output2Gain * output1Gain);
+            } else {
+                const double washoutGain = output3Gain * output2Gain * output1Gain * washoutScale;
+                matrixData.assign(refAlg, refDiff + transducer1State, washoutGain);
+                matrixData.assign(refAlg, refDiff + transducer2State, washoutGain);
+                matrixData.assign(refAlg, refDiff + washoutState, -washoutGain);
+            }
+        }
+    }
+    if (!hasDifferential(sMode)) {
+        return;
+    }
+
+    const auto addSelectedInput = [&matrixData, &inputLocs](index_t row, int mode, double value) {
+        if (mode == 1) {
+            matrixData.assignCheckCol(row, inputLocs[pssOmegaInLocation], value);
+        } else if (mode == 3) {
+            matrixData.assignCheckCol(row, inputLocs[pssElectricalPowerInLocation], value);
+        } else if (mode == 4) {
+            matrixData.assignCheckCol(row, inputLocs[pssPmechInLocation], value);
+        } else if (mode == 5) {
+            matrixData.assignCheckCol(row, inputLocs[pssVoltageInLocation], value);
+        }
+    };
+
+    matrixData.assign(refDiff + transducer1State,
+                      refDiff + transducer1State,
+                      ((-1.0 / T1) - stateData.cj));
+    addSelectedInput(refDiff + transducer1State, mode1, K1 / T1);
+    matrixData.assign(refDiff + transducer2State,
+                      refDiff + transducer2State,
+                      ((-1.0 / T2) - stateData.cj));
+    addSelectedInput(refDiff + transducer2State, mode2, K2 / T2);
+    matrixData.assign(refDiff + washoutState, refDiff + transducer1State, 1.0 / T4);
+    matrixData.assign(refDiff + washoutState, refDiff + transducer2State, 1.0 / T4);
+    matrixData.assign(refDiff + washoutState, refDiff + washoutState, ((-1.0 / T4) - stateData.cj));
+
+    matrixData.assign(refDiff + leadLag1State,
+                      refDiff + leadLag1State,
+                      ((-1.0 / T6) - stateData.cj));
+    if (T3 <= 0.0) {
+        matrixData.assign(refDiff + leadLag1State, refDiff + washoutState, 1.0 / T6);
+    } else {
+        const double gain = washoutScale / T6;
+        matrixData.assign(refDiff + leadLag1State, refDiff + transducer1State, gain);
+        matrixData.assign(refDiff + leadLag1State, refDiff + transducer2State, gain);
+        matrixData.assign(refDiff + leadLag1State, refDiff + washoutState, -gain);
+    }
+
+    matrixData.assign(refDiff + leadLag2State,
+                      refDiff + leadLag2State,
+                      ((-1.0 / T8) - stateData.cj));
+    matrixData.assign(refDiff + leadLag2State, refDiff + leadLag1State, ll1StateGain / T8);
+    if (T3 <= 0.0) {
+        matrixData.assign(refDiff + leadLag2State, refDiff + washoutState, output1Gain / T8);
+    } else {
+        const double gain = output1Gain * washoutScale / T8;
+        matrixData.assign(refDiff + leadLag2State, refDiff + transducer1State, gain);
+        matrixData.assign(refDiff + leadLag2State, refDiff + transducer2State, gain);
+        matrixData.assign(refDiff + leadLag2State, refDiff + washoutState, -gain);
+    }
+
+    matrixData.assign(refDiff + leadLag3State,
+                      refDiff + leadLag3State,
+                      ((-1.0 / T10) - stateData.cj));
+    matrixData.assign(refDiff + leadLag3State, refDiff + leadLag2State, ll2StateGain / T10);
+    matrixData.assign(refDiff + leadLag3State,
+                      refDiff + leadLag1State,
+                      output2Gain * ll1StateGain / T10);
+    if (T3 <= 0.0) {
+        matrixData.assign(refDiff + leadLag3State,
+                          refDiff + washoutState,
+                          output2Gain * output1Gain / T10);
+    } else {
+        const double gain = output2Gain * output1Gain * washoutScale / T10;
+        matrixData.assign(refDiff + leadLag3State, refDiff + transducer1State, gain);
+        matrixData.assign(refDiff + leadLag3State, refDiff + transducer2State, gain);
+        matrixData.assign(refDiff + leadLag3State, refDiff + washoutState, -gain);
+    }
+}
+
+void StabilizerST2CUT::timestep(CoreTime time, const IOdata& inputs, const SolverMode& /*sMode*/)
+{
+    derivative(inputs, emptyStateData, m_dstate_dt.data(), cLocalSolverMode);
+    const double timeStep = time - prevTime;
+    double* state = m_state.data() + 1;
+    const double* stateDerivative = m_dstate_dt.data() + 1;
+    for (index_t index = 0; index < 6; ++index) {
+        state[index] += timeStep * stateDerivative[index];
+    }
+    m_state[0] = output(inputs, state);
+    updateLimitFlags(inputs, state);
+    prevTime = time;
+}
+
+void StabilizerST2CUT::rootTest(const IOdata& inputs,
+                                const StateData& stateData,
+                                double roots[],
+                                const SolverMode& sMode)
+{
+    const auto locations = offsets.getLocations(stateData, sMode, this);
+    const index_t rootOffset = offsets.getRootOffset(sMode);
+    const double unlimited = unlimitedOutput(locations.diffStateLoc);
+    roots[rootOffset] = Lsmax - unlimited;
+    roots[rootOffset + 1] = unlimited - Lsmin;
+    roots[rootOffset + 2] = inputs[pssVoltageInLocation] - (initialVoltage + Vcl);
+    roots[rootOffset + 3] = (initialVoltage + Vcu) - inputs[pssVoltageInLocation];
+}
+
+void StabilizerST2CUT::rootTrigger(CoreTime /*time*/,
+                                   const IOdata& inputs,
+                                   const std::vector<int>& rootMask,
+                                   const SolverMode& sMode)
+{
+    const index_t rootOffset = offsets.getRootOffset(sMode);
+    if ((rootMask[rootOffset] == 0) && (rootMask[rootOffset + 1] == 0) &&
+        (rootMask[rootOffset + 2] == 0) && (rootMask[rootOffset + 3] == 0)) {
+        return;
+    }
+    if (updateLimitFlags(inputs, m_state.data() + 1)) {
+        alert(this, JAC_COUNT_CHANGE);
+    }
+}
+
+ChangeCode StabilizerST2CUT::rootCheck(const IOdata& inputs,
+                                       const StateData& /*stateData*/,
+                                       const SolverMode& /*sMode*/,
+                                       CheckLevel /*level*/)
+{
+    if (updateLimitFlags(inputs, m_state.data() + 1)) {
+        alert(this, JAC_COUNT_CHANGE);
+        return ChangeCode::JACOBIAN_CHANGE;
+    }
+    return ChangeCode::NO_CHANGE;
+}
+
+stringVec StabilizerST2CUT::localStateNames() const
+{
+    return {"vs", "l1", "l2", "wo", "ll1", "ll2", "ll3"};
+}
+
+index_t StabilizerST2CUT::findIndex(std::string_view field, const SolverMode& sMode) const
+{
+    if ((field == "vss") || (field == "vs")) {
+        return getOutputLoc(sMode, 0);
+    }
+    return kInvalidLocation;
+}
+
+void StabilizerST2CUT::set(std::string_view param, std::string_view val)
+{
+    Stabilizer::set(param, val);
+}
+
+void StabilizerST2CUT::set(std::string_view param, double val, units::unit unitType)
+{
+    const auto setFinite = [val](const char* parameterName) {
+        if (!std::isfinite(val)) {
+            throw InvalidParameterValue(std::string("ST2CUT ") + parameterName + " must be finite");
+        }
+    };
+    if ((param == "mode") || (param == "mode1")) {
+        if (!std::isfinite(val) || (std::floor(val) != val) ||
+            !supportedMode(static_cast<int>(val))) {
+            throw InvalidParameterValue("ST2CUT MODE is unsupported");
+        }
+        mode1 = static_cast<int>(val);
+    } else if ((param == "mode2") || (param == "mode_2")) {
+        if (!std::isfinite(val) || (std::floor(val) != val) ||
+            !supportedMode(static_cast<int>(val))) {
+            throw InvalidParameterValue("ST2CUT MODE2 is unsupported");
+        }
+        mode2 = static_cast<int>(val);
+    } else if ((param == "busr") || (param == "busr1")) {
+        if (!std::isfinite(val) || (std::floor(val) != val) || (val != 0.0)) {
+            throw InvalidParameterValue("ST2CUT remote BUSR is unsupported");
+        }
+        remoteBus1 = 0;
+    } else if (param == "busr2") {
+        if (!std::isfinite(val) || (std::floor(val) != val) || (val != 0.0)) {
+            throw InvalidParameterValue("ST2CUT remote BUSR2 is unsupported");
+        }
+        remoteBus2 = 0;
+    } else if (param == "k1") {
+        setFinite("K1");
+        K1 = val;
+    } else if (param == "k2") {
+        setFinite("K2");
+        K2 = val;
+    } else if (param == "t1") {
+        setFinite("T1");
+        if (val <= 0.0) {
+            throw InvalidParameterValue("ST2CUT T1 must be positive");
+        }
+        T1 = val;
+    } else if (param == "t2") {
+        setFinite("T2");
+        if (val <= 0.0) {
+            throw InvalidParameterValue("ST2CUT T2 must be positive");
+        }
+        T2 = val;
+    } else if (param == "t3") {
+        setFinite("T3");
+        if (val < 0.0) {
+            throw InvalidParameterValue("ST2CUT T3 must be nonnegative");
+        }
+        T3 = val;
+    } else if (param == "t4") {
+        setFinite("T4");
+        if (val <= 0.0) {
+            throw InvalidParameterValue("ST2CUT T4 must be positive");
+        }
+        T4 = val;
+    } else if (param == "t5") {
+        setFinite("T5");
+        if (val < 0.0) {
+            throw InvalidParameterValue("ST2CUT T5 must be nonnegative");
+        }
+        T5 = val;
+    } else if (param == "t6") {
+        setFinite("T6");
+        if (val <= 0.0) {
+            throw InvalidParameterValue("ST2CUT T6 must be positive");
+        }
+        T6 = val;
+    } else if (param == "t7") {
+        setFinite("T7");
+        if (val < 0.0) {
+            throw InvalidParameterValue("ST2CUT T7 must be nonnegative");
+        }
+        T7 = val;
+    } else if (param == "t8") {
+        setFinite("T8");
+        if (val <= 0.0) {
+            throw InvalidParameterValue("ST2CUT T8 must be positive");
+        }
+        T8 = val;
+    } else if (param == "t9") {
+        setFinite("T9");
+        if (val < 0.0) {
+            throw InvalidParameterValue("ST2CUT T9 must be nonnegative");
+        }
+        T9 = val;
+    } else if (param == "t10") {
+        setFinite("T10");
+        if (val <= 0.0) {
+            throw InvalidParameterValue("ST2CUT T10 must be positive");
+        }
+        T10 = val;
+    } else if ((param == "lsmax") || (param == "vmax")) {
+        setFinite("LSMAX");
+        if (val < Lsmin) {
+            throw InvalidParameterValue("ST2CUT LSMAX must not be less than LSMIN");
+        }
+        Lsmax = val;
+    } else if ((param == "lsmin") || (param == "vmin")) {
+        setFinite("LSMIN");
+        if (val > Lsmax) {
+            throw InvalidParameterValue("ST2CUT LSMIN must not exceed LSMAX");
+        }
+        Lsmin = val;
+    } else if (param == "vcu") {
+        setFinite("VCU");
+        const double mappedValue = (val == 0.0) ? 999.0 : val;
+        if (mappedValue < Vcl) {
+            throw InvalidParameterValue("ST2CUT VCU must not be less than VCL");
+        }
+        Vcu = mappedValue;
+    } else if (param == "vcl") {
+        setFinite("VCL");
+        const double mappedValue = (val == 0.0) ? -999.0 : val;
+        if (mappedValue > Vcu) {
+            throw InvalidParameterValue("ST2CUT VCL must not exceed VCU");
+        }
+        Vcl = mappedValue;
+    } else {
+        Stabilizer::set(param, val, unitType);
+    }
+}
+
+double StabilizerST2CUT::get(std::string_view param, units::unit unitType) const
+{
+    if ((param == "mode") || (param == "mode1")) {
+        return mode1;
+    }
+    if (param == "mode2") {
+        return mode2;
+    }
+    if ((param == "busr") || (param == "busr1")) {
+        return remoteBus1;
+    }
+    if (param == "busr2") {
+        return remoteBus2;
+    }
+    if (param == "k1") {
+        return K1;
+    }
+    if (param == "k2") {
+        return K2;
+    }
+    if (param == "t1") {
+        return T1;
+    }
+    if (param == "t2") {
+        return T2;
+    }
+    if (param == "t3") {
+        return T3;
+    }
+    if (param == "t4") {
+        return T4;
+    }
+    if (param == "t5") {
+        return T5;
+    }
+    if (param == "t6") {
+        return T6;
+    }
+    if (param == "t7") {
+        return T7;
+    }
+    if (param == "t8") {
+        return T8;
+    }
+    if (param == "t9") {
+        return T9;
+    }
+    if (param == "t10") {
+        return T10;
+    }
+    if ((param == "lsmax") || (param == "vmax")) {
+        return Lsmax;
+    }
+    if ((param == "lsmin") || (param == "vmin")) {
+        return Lsmin;
+    }
+    if (param == "vcu") {
+        return Vcu;
+    }
+    if (param == "vcl") {
+        return Vcl;
+    }
+    return Stabilizer::get(param, unitType);
+}
+}  // namespace griddyn::stabilizers
diff --git a/src/griddyn/stabilizers/StabilizerST2CUT.h b/src/griddyn/stabilizers/StabilizerST2CUT.h
new file mode 100644
index 000000000..57c1b123f
--- /dev/null
+++ b/src/griddyn/stabilizers/StabilizerST2CUT.h
@@ -0,0 +1,179 @@
+/*
+ * Copyright (c) 2014-2026, Lawrence Livermore National Security
+ * See the top-level NOTICE for additional details. All rights reserved.
+ * SPDX-License-Identifier: BSD-3-Clause
+ */
+
+#pragma once
+
+#include "../Stabilizer.h"
+#include <string>
+#include <vector>
+
+namespace griddyn::stabilizers {
+/**
+ * @brief PSS/E ST2CUT dual-input power-system stabilizer.
+ *
+ * This is the frozen ANDES v2.0.0 realization in
+ * `andes/models/pss/st2cut.py`, using the PSS/E DYR field order declared in
+ * `andes/io/psse-dyr.yaml`. It implements two transducer lags, their sum, a
+ * washout (or lag when @f$T_3=0@f$), three lead-lag stages, the
+ * @f$[LSMIN,LSMAX]@f$ output limiter, and the initialized-voltage-relative
+ * @f$[VCL,VCU]@f$ output gate.
+ *
+ * Let @f$u_1,u_2@f$ be the selected measurement signals; @f$x_1,x_2@f$ the
+ * transducer states; @f$x_w@f$ the washout/lag state; and
+ * @f$x_5,x_7,x_9@f$ the three lead-lag lag states. The local differential
+ * equations are
+ *
+ * @f[
+ * T_1\dot{x}_1=K_1u_1-x_1,\qquad T_2\dot{x}_2=K_2u_2-x_2,\qquad
+ * T_4\dot{x}_w=x_1+x_2-x_w.
+ * @f]
+ *
+ * The washout-or-lag output is
+ *
+ * @f[
+ * w=\begin{cases}
+ * \dfrac{T_3}{T_4}(x_1+x_2-x_w), & T_3>0,\\
+ * x_w, & T_3=0,
+ * \end{cases}
+ * @f]
+ *
+ * and each lead-lag block @f$(T_a,T_b,x,y)@f$ obeys
+ *
+ * @f[
+ * T_b\dot{x}=u-x,\qquad y=x+\dfrac{T_a}{T_b}(u-x).
+ * @f]
+ *
+ * The three blocks are cascaded as
+ * @f$(w,T_5,T_6,x_5,y_5)@f$, @f$(y_5,T_7,T_8,x_7,y_7)@f$, and
+ * @f$(y_7,T_9,T_{10},x_9,y_9)@f$. Finally, with initialized terminal-voltage
+ * reference @f$V_{t0}@f$,
+ *
+ * @f[
+ * V_{SS}=\begin{cases}
+ * \operatorname{limit}_{[LSMIN,LSMAX]}(y_9), &
+ * V_{t0}+VCL\leq V_t\leq V_{t0}+VCU,\\
+ * 0, & \text{otherwise}.
+ * \end{cases}
+ * @f]
+ *
+ * Initialization sets @f$x_1=K_1u_1@f$, @f$x_2=K_2u_2@f$,
+ * @f$x_w=x_1+x_2@f$, and the three lead-lag states and @f$V_{SS}@f$ to zero,
+ * so the washout path contributes no steady-state stabilizing signal.
+ *
+ * @par References
+ * - ANDES v2.0.0, `andes/models/pss/st2cut.py` (equations and initialization).
+ * - ANDES v2.0.0, `andes/io/psse-dyr.yaml` (DYR schema and field order).
+ * - Siemens PTI, <em>PSS/E Model Library</em>, ST2CUT model (published PSS/E
+ *   block-model source).
+ *
+ * The DYR adapter maps the exact fields
+ * `BUS, ID, MODE, BUSR, MODE2, BUSR2, K1, K2, T1, T2, T3, T4, T5, T6, T7,
+ * T8, T9, T10, LSMAX, LSMIN, VCU, VCL`. Following ANDES, a zero `VCU` or
+ * `VCL` is expanded to `+999` or `-999`, respectively, to disable that side
+ * of the voltage gate.
+ *
+ * The GridDyn PSS signal contract exposes rotor speed, terminal voltage,
+ * mechanical power, and electrical power.  It therefore supports local
+ * MODE/MODE2 values 0 (disabled), 1 (rotor speed), 3 (electrical power),
+ * 4 (accelerating power), and 5 (terminal voltage).  ANDES modes 2
+ * (BusFreq) and 6 (voltage derivative), and nonzero remote BUSR/BUSR2,
+ * require cross-bus measurement routing that GridDyn does not yet provide;
+ * they are rejected at initialization rather than silently approximated.
+ * In particular, MODE 1 is @f$\omega-1@f$, MODE 3 is electrical power,
+ * MODE 4 is @f$P_m-P_{m0}@f$, and MODE 5 is terminal voltage. MODE 0 supplies
+ * zero to its transducer.
+ *
+ * The implementation requires positive @f$T_1,T_2,T_4,T_6,T_8,T_{10}@f$.
+ * Nonnegative @f$T_3,T_5,T_7,T_9@f$ are accepted; @f$T_3=0@f$ selects the
+ * documented lag form above. Other zero-denominator block bypasses are not
+ * implemented and are rejected during parameter validation.
+ */
+class StabilizerST2CUT: public Stabilizer {
+  public:
+    enum ST2CUTFlags {
+        OUTPUT_LIMITED = OBJECT_FLAG5,
+        OUTPUT_LIMIT_HIGH = OBJECT_FLAG6,
+        VOLTAGE_GATED = OBJECT_FLAG7,
+    };
+
+  protected:
+    int mode1 = 1;
+    int remoteBus1 = 0;
+    int mode2 = 0;
+    int remoteBus2 = 0;
+    double K1 = 1.0;
+    double K2 = 1.0;
+    double T1 = 1.0;
+    double T2 = 1.0;
+    double T3 = 1.0;
+    double T4 = 0.2;
+    double T5 = 1.0;
+    double T6 = 0.5;
+    double T7 = 1.0;
+    double T8 = 1.0;
+    double T9 = 1.0;
+    double T10 = 0.2;
+    double Lsmax = 0.3;
+    double Lsmin = -0.3;
+    double Vcu = 999.0;
+    double Vcl = -999.0;
+    double initialVoltage = 1.0;
+    double initialPmech = 0.0;
+
+  public:
+    explicit StabilizerST2CUT(const std::string& objName = "pssST2CUT_#");
+    CoreObject* clone(CoreObject* obj = nullptr) const override;
+
+    void dynObjectInitializeA(CoreTime time0, std::uint32_t flags) override;
+    void dynObjectInitializeB(const IOdata& inputs,
+                              const IOdata& desiredOutput,
+                              IOdata& fieldSet) override;
+
+    void set(std::string_view param, std::string_view val) override;
+    void set(std::string_view param, double val, units::unit unitType = units::defunit) override;
+    double get(std::string_view param, units::unit unitType = units::defunit) const override;
+
+    stringVec localStateNames() const override;
+    index_t findIndex(std::string_view field, const SolverMode& sMode) const override;
+    void residual(const IOdata& inputs,
+                  const StateData& stateData,
+                  double resid[],
+                  const SolverMode& sMode) override;
+    void derivative(const IOdata& inputs,
+                    const StateData& stateData,
+                    double deriv[],
+                    const SolverMode& sMode) override;
+    void jacobianElements(const IOdata& inputs,
+                          const StateData& stateData,
+                          MatrixData<double>& matrixData,
+                          const IOlocs& inputLocs,
+                          const SolverMode& sMode) override;
+    void timestep(CoreTime time, const IOdata& inputs, const SolverMode& sMode) override;
+    void rootTest(const IOdata& inputs,
+                  const StateData& stateData,
+                  double roots[],
+                  const SolverMode& sMode) override;
+    void rootTrigger(CoreTime time,
+                     const IOdata& inputs,
+                     const std::vector<int>& rootMask,
+                     const SolverMode& sMode) override;
+    ChangeCode rootCheck(const IOdata& inputs,
+                         const StateData& stateData,
+                         const SolverMode& sMode,
+                         CheckLevel level) override;
+
+  private:
+    double selectedInput(const IOdata& inputs, int mode) const;
+    double washoutOutput(const double state[]) const;
+    static double leadLagOutput(double input, double state, double leadTime, double lagTime);
+    double unlimitedOutput(const double state[]) const;
+    bool voltageEnabled(const IOdata& inputs) const;
+    double output(const IOdata& inputs, const double state[]) const;
+    int outputLimitStatus(const double state[]) const;
+    bool updateLimitFlags(const IOdata& inputs, const double state[]);
+    static bool supportedMode(int mode);
+};
+}  // namespace griddyn::stabilizers
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 1d0a82c77..9f186392a 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -57,6 +57,7 @@ set(generator_component_gtest_sources
     gtestHelperFunctions.cpp
     componentTests/testBlocks.cpp
     componentTests/testGovernors.cpp
+    componentTests/testStabilizers.cpp
     componentTests/testGenerators.cpp
     componentTests/testGenModels.cpp
     componentTests/testExciters.cpp
diff --git a/test/andesTests/testAndesDyrReader.cpp b/test/andesTests/testAndesDyrReader.cpp
index e87d88859..fd7d3ce65 100644
--- a/test/andesTests/testAndesDyrReader.cpp
+++ b/test/andesTests/testAndesDyrReader.cpp
@@ -17,6 +17,7 @@
 #include "griddyn/genmodels/GenModelGENROU.h"
 #include "griddyn/governors/GovernorIeeeG1.h"
 #include "griddyn/governors/GovernorTgov1.h"
+#include "griddyn/stabilizers/StabilizerST2CUT.h"
 #include <algorithm>
 #include <array>
 #include <cmath>
@@ -27,6 +28,7 @@
 #include <nlohmann/json.hpp>
 #include <string>
 #include <string_view>
+#include <tuple>
 #include <utility>
 #include <vector>
 
@@ -75,7 +77,7 @@ std::vector<double> runLoadStepCase(const std::vector<std::string_view>& dyrFile
         if (generator == nullptr) {
             continue;
         }
-        for (const auto controllerName : {"governor", "exciter"}) {
+        for (const auto controllerName : {"governor", "exciter", "pss"}) {
             auto* controller =
                 dynamic_cast<griddyn::GridSubModel*>(generator->find(controllerName));
             if (controller != nullptr) {
@@ -361,6 +363,34 @@ TEST(AndesDyrReaderTests, MapsExst1ParametersAndCouplesToGenrou)
     EXPECT_EQ(runJacobianCheck(simulation, griddyn::cDaeSolverMode, false), 0);
 }
 
+TEST(AndesDyrReaderTests, MapsSt2cutParametersAndCouplesToExciters)
+{
+    auto simulation = std::make_unique<griddyn::GridDynSimulation>();
+    griddyn::loadFile(simulation.get(), makeAndesTestPath("ieee14.raw"));
+    griddyn::loadFile(simulation.get(), makeAndesTestPath("ieee14_genrou.dyr"));
+    griddyn::loadFile(simulation.get(), makeAndesTestPath("ieee14_esst3a.dyr"));
+    griddyn::loadFile(simulation.get(), makeAndesTestPath("ieee14_exst1.dyr"));
+    griddyn::loadFile(simulation.get(), makeAndesTestPath("ieee14_st2cut.dyr"));
+
+    for (const auto& [busId, expectedK1, expectedLsmax] :
+         {std::tuple{1, 1.2, 0.05}, std::tuple{2, 1.1, 0.06}}) {
+        auto* bus = dynamic_cast<griddyn::GridBus*>(simulation->findByUserID("bus", busId));
+        ASSERT_NE(bus, nullptr);
+        auto* generator = dynamic_cast<griddyn::DynamicGenerator*>(bus->getGen(0));
+        ASSERT_NE(generator, nullptr);
+        auto* stabilizer =
+            dynamic_cast<griddyn::stabilizers::StabilizerST2CUT*>(generator->find("pss"));
+        ASSERT_NE(stabilizer, nullptr);
+        EXPECT_DOUBLE_EQ(stabilizer->get("mode"), 1.0);
+        EXPECT_DOUBLE_EQ(stabilizer->get("mode2"), 0.0);
+        EXPECT_DOUBLE_EQ(stabilizer->get("k1"), expectedK1);
+        EXPECT_DOUBLE_EQ(stabilizer->get("lsmax"), expectedLsmax);
+    }
+    EXPECT_EQ(simulation->dynInitialize(), 0);
+    EXPECT_EQ(runResidualCheck(simulation, griddyn::cDaeSolverMode, false), 0);
+    EXPECT_EQ(runJacobianCheck(simulation, griddyn::cDaeSolverMode, false), 0);
+}
+
 TEST(AndesDynamicTests, Tgov1RespondsToLoadStep)
 {
     const auto finalState = runLoadStepCase({"ieee14_tgov1.dyr"});
@@ -385,6 +415,13 @@ TEST(AndesDynamicTests, Exst1RespondsToLoadStep)
     EXPECT_FALSE(finalState.empty());
 }
 
+TEST(AndesDynamicTests, St2cutRespondsToLoadStep)
+{
+    const auto finalState =
+        runLoadStepCase({"ieee14_esst3a.dyr", "ieee14_exst1.dyr", "ieee14_st2cut.dyr"});
+    EXPECT_FALSE(finalState.empty());
+}
+
 TEST(AndesDynamicTests, CombinedControllersRespondToLoadStep)
 {
     const auto finalState = runLoadStepCase(
diff --git a/test/componentTests/testStabilizers.cpp b/test/componentTests/testStabilizers.cpp
new file mode 100644
index 000000000..94979af5b
--- /dev/null
+++ b/test/componentTests/testStabilizers.cpp
@@ -0,0 +1,150 @@
+/*
+ * Copyright (c) 2014-2026, Lawrence Livermore National Security
+ * See the top-level NOTICE for additional details. All rights reserved.
+ * SPDX-License-Identifier: BSD-3-Clause
+ */
+
+#include "../gtestHelper.h"
+#include "core/ObjectFactory.hpp"
+#include "griddyn/Stabilizer.h"
+#include "griddyn/generators/DynamicGenerator.h"
+#include "griddyn/stabilizers/StabilizerST2CUT.h"
+#include <gtest/gtest.h>
+#include <memory>
+#include <string>
+#include <vector>
+
+using namespace griddyn;
+
+class StabilizerTests: public GridDynSimulationTestFixture, public ::testing::Test {};
+
+namespace {
+void configureSt2cut(stabilizers::StabilizerST2CUT& stabilizer)
+{
+    stabilizer.set("mode", 1.0);
+    stabilizer.set("mode2", 0.0);
+    stabilizer.set("k1", 1.2);
+    stabilizer.set("k2", 0.8);
+    stabilizer.set("t1", 0.1);
+    stabilizer.set("t2", 0.2);
+    stabilizer.set("t3", 3.0);
+    stabilizer.set("t4", 1.0);
+    stabilizer.set("t5", 0.5);
+    stabilizer.set("t6", 0.1);
+    stabilizer.set("t7", 0.3);
+    stabilizer.set("t8", 0.15);
+    stabilizer.set("t9", 0.2);
+    stabilizer.set("t10", 0.1);
+    stabilizer.set("lsmax", 0.05);
+    stabilizer.set("lsmin", -0.05);
+    stabilizer.set("vcu", 0.1);
+    stabilizer.set("vcl", -0.1);
+}
+}  // namespace
+
+TEST(StabilizerModelTests, St2cutMatchesAndesInitializationAndPerturbedEquations)
+{
+    stabilizers::StabilizerST2CUT stabilizer;
+    configureSt2cut(stabilizer);
+    stabilizer.dynInitializeA(0.0, 0);
+    IOdata inputs{1.0, 1.0, 0.8, 0.7};
+    IOdata fieldSet;
+    stabilizer.dynInitializeB(inputs, {0.0}, fieldSet);
+
+    const auto& initialized = stabilizer.getStates();
+    ASSERT_EQ(initialized.size(), 7U);
+    for (const auto value : initialized) {
+        EXPECT_NEAR(value, 0.0, 1e-14);
+    }
+
+    // State order is [VSS, L1_y, L2_y, WO_x, LL1_x, LL2_x, LL3_x].
+    // These are direct evaluations of frozen ANDES ST2CUT equations.
+    std::vector<double> state{0.0, 0.02, -0.01, 0.005, 0.004, 0.003, 0.002};
+    std::vector<double> stateDerivative(state.size(), 0.0);
+    stabilizer.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode);
+    inputs[pssOmegaInLocation] = 1.01;
+
+    std::vector<double> derivative(state.size(), 0.0);
+    stabilizer.derivative(inputs, emptyStateData, derivative.data(), cLocalSolverMode);
+    EXPECT_NEAR(derivative[1], -0.08, 1e-14);
+    EXPECT_NEAR(derivative[2], 0.05, 1e-14);
+    EXPECT_NEAR(derivative[3], 0.005, 1e-14);
+    EXPECT_NEAR(derivative[4], 0.11, 1e-14);
+    EXPECT_NEAR(derivative[5], 0.37333333333333335, 1e-14);
+    EXPECT_NEAR(derivative[6], 1.13, 1e-14);
+
+    std::vector<double> residual(state.size(), 0.0);
+    stabilizer.residual(inputs, emptyStateData, residual.data(), cLocalSolverMode);
+    EXPECT_NEAR(residual[0], 0.05, 1e-14);
+    for (std::size_t index = 1; index < residual.size(); ++index) {
+        EXPECT_NEAR(residual[index], derivative[index], 1e-14) << index;
+    }
+}
+
+TEST(StabilizerModelTests, St2cutOutputLimitsVoltageGatingAndRejectsUnsupportedSignals)
+{
+    stabilizers::StabilizerST2CUT stabilizer;
+    configureSt2cut(stabilizer);
+    stabilizer.dynInitializeA(0.0, 0);
+    stabilizer.setRootOffset(0, cLocalSolverMode);
+    IOdata inputs{1.0, 1.0, 0.8, 0.7};
+    IOdata fieldSet;
+    stabilizer.dynInitializeB(inputs, {0.0}, fieldSet);
+
+    std::vector<double> state{0.0, 0.02, -0.01, 0.005, 0.004, 0.003, 0.002};
+    std::vector<double> stateDerivative(state.size(), 0.0);
+    stabilizer.setState(0.0, state.data(), stateDerivative.data(), cLocalSolverMode);
+    std::vector<double> residual(state.size(), 0.0);
+    stabilizer.residual(inputs, emptyStateData, residual.data(), cLocalSolverMode);
+    EXPECT_NEAR(residual[0], 0.05, 1e-14);
+
+    inputs[pssVoltageInLocation] = 1.2;
+    stabilizer.residual(inputs, emptyStateData, residual.data(), cLocalSolverMode);
+    EXPECT_NEAR(residual[0], 0.0, 1e-14);
+
+    std::vector<double> roots(4, 0.0);
+    stabilizer.rootTest(inputs, emptyStateData, roots.data(), cLocalSolverMode);
+    EXPECT_LT(roots[0], 0.0);
+    EXPECT_LT(roots[3], 0.0);
+    stabilizer.rootTrigger(0.0, inputs, {1, 0, 0, 0}, cLocalSolverMode);
+    EXPECT_TRUE(stabilizer.checkFlag(stabilizers::StabilizerST2CUT::OUTPUT_LIMITED));
+    EXPECT_TRUE(stabilizer.checkFlag(stabilizers::StabilizerST2CUT::OUTPUT_LIMIT_HIGH));
+    stabilizer.rootTrigger(0.0, inputs, {0, 0, 0, 1}, cLocalSolverMode);
+    EXPECT_TRUE(stabilizer.checkFlag(stabilizers::StabilizerST2CUT::VOLTAGE_GATED));
+
+    EXPECT_ANY_THROW(stabilizer.set("mode", 2.0));
+    EXPECT_ANY_THROW(stabilizer.set("mode2", 6.0));
+    EXPECT_ANY_THROW(stabilizer.set("busr", 9.0));
+}
+
+TEST(StabilizerModelTests, St2cutFactoryCloneAndParameterValidation)
+{
+    auto factory = CoreObjectFactory::instance();
+    std::unique_ptr<CoreObject> object(factory->createObject("pss", "st2cut"));
+    auto* stabilizer = dynamic_cast<stabilizers::StabilizerST2CUT*>(object.get());
+    ASSERT_NE(stabilizer, nullptr);
+    stabilizer->set("k1", 2.5);
+    std::unique_ptr<CoreObject> clonedObject(stabilizer->clone());
+    auto* clone = dynamic_cast<stabilizers::StabilizerST2CUT*>(clonedObject.get());
+    ASSERT_NE(clone, nullptr);
+    EXPECT_DOUBLE_EQ(clone->get("k1"), 2.5);
+
+    EXPECT_ANY_THROW(stabilizer->set("t4", 0.0));
+    EXPECT_ANY_THROW(stabilizer->set("lsmin", 0.4));
+}
+
+TEST_F(StabilizerTests, St2cutAnalyticJacobianMatchesFiniteDifferencesWhenAttached)
+{
+    gds = readSimXMLFile(std::string(GRIDDYN_TEST_DIRECTORY "/genmodel_tests/test_model1.xml"));
+    auto* generator = dynamic_cast<DynamicGenerator*>(gds->getGen(0));
+    ASSERT_NE(generator, nullptr);
+    auto* stabilizer = new stabilizers::StabilizerST2CUT();
+    configureSt2cut(*stabilizer);
+    stabilizer->set("lsmax", 2.0);
+    stabilizer->set("lsmin", -2.0);
+    generator->add(stabilizer);
+
+    ASSERT_EQ(gds->dynInitialize(), 0);
+    EXPECT_EQ(runResidualCheck(gds, cDaeSolverMode, false), 0);
+    EXPECT_EQ(runJacobianCheck(gds, cDaeSolverMode, false), 0);
+}
diff --git a/test/test_files/andes_tests/ieee14_st2cut.dyr b/test/test_files/andes_tests/ieee14_st2cut.dyr
new file mode 100644
index 000000000..0d6e4e75f
--- /dev/null
+++ b/test/test_files/andes_tests/ieee14_st2cut.dyr
@@ -0,0 +1,4 @@
+      1 'ST2CUT' 1 1 0 0 0 1.20 0.80 0.10 0.20 3.00 1.00 0.50 0.10
+          0.30 0.15 0.20 0.10 0.05 -0.05 0.0 0.0 /
+      2 'ST2CUT' 1 1 0 0 0 1.10 0.70 0.12 0.22 2.50 0.90 0.40 0.12
+          0.25 0.16 0.18 0.11 0.06 -0.06 0.0 0.0 /