add and test SCRX model

hexaeder · hexaeder · commit 4e0f43f42118 · 2025-09-27T16:25:33.000+02:00
diff --git a/src/Library/Library.jl b/src/Library/Library.jl
@@ -174,6 +174,9 @@ include("OpenIPSL/Loads/PSSE_Load.jl")
 export PSSE_IEEET1, PSSE_ESST4B, PSSE_EXST1, PSSE_ESST1A
 include("OpenIPSL/Controls/PSSE_Excitation.jl")
 
+export PSSE_SCRX
+include("OpenIPSL/Controls/PSSE_SCRX.jl")
+
 export PSSE_IEEEG1
 include("OpenIPSL/Controls/PSSE_TurbineGovernors.jl")
 
diff --git a/src/Library/OpenIPSL/Controls/PSSE_SCRX.jl b/src/Library/OpenIPSL/Controls/PSSE_SCRX.jl
@@ -0,0 +1,100 @@
+# PSSE Static Excitation System Components
+#
+# Original work Copyright (c) 2016-2022 Luigi Vanfretti, ALSETLab, and contributors
+# Original work licensed under BSD 3-Clause License
+# Original source: https://github.com/OpenIPSL/OpenIPSL
+#
+# This Julia/PowerDynamics port maintains the same mathematical formulation
+# while adapting to PowerDynamics/ModelingToolkit framework conventions.
+
+@mtkmodel PSSE_SCRX begin
+    @structural_parameters begin
+        # Optional input enabling flags (following naming convention from other exciter models)
+        vothsg_input=false # Other signal input (typically zero if disabled)
+        vuel_input=false   # Under-excitation limiter input (typically zero if disabled)
+        voel_input=false   # Over-excitation limiter input (typically zero if disabled)
+        C_SWITCH=false     # Feeding selection: false for bus fed, true for solid fed
+    end
+
+    @parameters begin
+        # Excitation system parameters with OpenIPSL defaults
+        T_AT_B=0.1, [description="Ratio between regulator numerator (lead) and denominator (lag) time constants [pu]"]
+        T_B=1, [description="Regulator denominator (lag) time constant [s]"]
+        K=100, [description="Excitation power source output gain [pu]"]
+        T_E=0.005, [description="Excitation power source output time constant [s]"]
+        E_MIN=-10, [description="Minimum exciter output [pu]"]
+        E_MAX=10, [description="Maximum exciter output [pu]"]
+        r_cr_fd=10, [description="Ratio between crowbar circuit resistance and field circuit resistance [pu]"]
+
+        # Initialization parameters (determined during initialization)
+        V_REF, [guess=1, description="Voltage reference setpoint [pu]"]
+    end
+
+    @components begin
+        # Required inputs/outputs
+        ECOMP_in = RealInput()   # Terminal voltage measurement input
+        XADIFD_in = RealInput()  # Machine field current input
+        EFD_out = RealOutput()   # Field voltage output to generator
+
+        # Optional inputs (created conditionally based on structural parameters)
+        if vothsg_input
+            VOTHSG_in = RealInput()  # Other signal input
+        end
+        if vuel_input
+            VUEL_in = RealInput()    # Under-excitation limiter input
+        end
+        if voel_input
+            VOEL_in = RealInput()    # Over-excitation limiter input
+        end
+
+        # Building block components
+        leadlag = LeadLag(K=1, T1=T_AT_B*T_B, T2=T_B, guess=1)
+        amplifier = SimpleLagLim(K=K, T=T_E, outMin=E_MIN, outMax=E_MAX, guess=1)
+    end
+
+    @variables begin
+        # Signal processing variables
+        voltage_error(t), [description="Voltage error signal [pu]"]
+        sum_signal(t), [description="Combined signal before leadlag [pu]"]
+
+        # Amplifier and switching variables
+        amplifier_output(t), [description="Amplifier output before switching [pu]"]
+        switch_output(t), [description="Switched output (bus fed or solid fed) [pu]"]
+
+        # NegCurLogic variables
+        crowbar_voltage(t), [description="Crowbar circuit voltage [pu]"]
+        field_current(t), [description="Field current input [pu]"]
+        EFD(t), [description="Final exciter output [pu]"]
+    end
+
+    @equations begin
+        # Signal combination following OpenIPSL Add3 and DiffV1 logic
+        voltage_error ~ V_REF - ECOMP_in.u
+        sum_signal ~ voltage_error + (vothsg_input ? VOTHSG_IN.u : 0) + (vuel_input ? VUEL_IN.u : 0) - (voel_input ? VOEL_IN.u : 0)
+
+        # LeadLag compensator
+        leadlag.in ~ sum_signal
+
+        # Amplifier with limits
+        amplifier.in ~ leadlag.out
+        amplifier_output ~ amplifier.out
+
+        if C_SWITCH
+            # Solid fed: use amplifier output directly
+            switch_output ~ amplifier_output
+        else
+            # Bus fed: multiply by terminal voltage
+            switch_output ~ amplifier_output * ECOMP_in.u
+        end
+
+        # Field current input
+        field_current ~ XADIFD_in.u
+
+        # Integrated NegCurLogic functionality
+        crowbar_voltage ~ ifelse(abs(r_cr_fd) < eps(), 0.0, -r_cr_fd * field_current)
+        EFD ~ ifelse(field_current < 0, crowbar_voltage, switch_output)
+
+        # Output connection
+        EFD_out.u ~ EFD
+    end
+end
diff --git a/src/Library/building_blocks.jl b/src/Library/building_blocks.jl
@@ -170,6 +170,7 @@ function _generate_limint_callbacks(cf::NetworkDynamics.ComponentModel, namespac
         end
     end
 
+    # TODO: merge both discrete conditions and move condition below function for performance
     discrete_condition = ComponentCondition([out, min, max], []) do u, p, t
         # account for nummerical innaccuracies at the boudaries
         u[out] < u[min] - 1e-10 || u[out] > u[max] + 1e-10
@@ -199,10 +200,29 @@ function _generate_limint_callbacks(cf::NetworkDynamics.ComponentModel, namespac
             error("Sanity check was wrongfully triggered!")
         end
     end
+    discrete_unsat_condition = ComponentCondition([forcing],[satmin, satmax]) do u, p, t
+        insatmin = !iszero(p[satmin])
+        insatmax = !iszero(p[satmax])
+        insatmin && u[forcing] > 0 || insatmax && u[forcing] < 0
+    end
+    discrete_unsat_affect = ComponentAffect([],[satmin, satmax]) do u, p, ctx
+        insatmin = !iszero(p[satmin])
+        insatmax = !iszero(p[satmax])
+        if insatmin
+            println("$namespace: _/ left lower saturation at $(round(ctx.t, digits=4))s (triggered by discrete cb)")
+            p[satmin] = 0.0
+        elseif insatmax
+            println("$namespace: ⎺\\ left upper saturation at $(round(ctx.t, digits=4))s (triggered by discrete cb)")
+            p[satmax] = 0.0
+        else
+            error("Sanity check was wrongfully triggered!")
+        end
+    end
 
     (
         VectorContinuousComponentCallback(condition, upcrossing_affect, 3; affect_neg! = downcrossing_affect),
-        DiscreteComponentCallback(discrete_condition, discrete_affect)
+        DiscreteComponentCallback(discrete_condition, discrete_affect),
+        DiscreteComponentCallback(discrete_unsat_condition, discrete_unsat_affect)
     )
 end
 function _SatLim_condition(_out, u, p, _)
diff --git a/test/OpenIPSL_test/ModelTransferRules.md b/test/OpenIPSL_test/ModelTransferRules.md
@@ -212,7 +212,7 @@ The following OpenIPSL models use the SMIB test harness and are candidates for P
   - *Dependencies: Requires GENROE machine model*
 - [ ] **SEXS** - Simplified excitation system
   - *Dependencies: Requires GENROU machine model*
-- [ ] **SCRX** - Static exciter
+- [X] **SCRX** - Static exciter
   - *Dependencies: Requires GENROU machine model*
 - [X] **EXST1** - Static excitation system
   - *Dependencies: Requires GENROE machine model*
diff --git a/test/OpenIPSL_test/PSSE_SCRX_test.jl b/test/OpenIPSL_test/PSSE_SCRX_test.jl
@@ -0,0 +1,141 @@
+using PowerDynamics
+PowerDynamics.load_pdtesting()
+using Main.PowerDynamicsTesting
+
+using PowerDynamics.Library
+using ModelingToolkit
+using OrdinaryDiffEqRosenbrock
+using OrdinaryDiffEqNonlinearSolve
+
+using CSV
+using DataFrames
+using CairoMakie
+
+# Load reference data from SCRX simulation
+ref = CSV.read(
+    joinpath(pkgdir(PowerDynamics),"test","OpenIPSL_test","SCRX","modelica_results.csv.gz"),
+    DataFrame;
+    drop=(i,name) -> contains(string(name), "nrows="),
+    silencewarnings=true
+)
+
+# SCRX test system parameters from OpenIPSL SCRX test case
+BUS = let
+    # Power flow from OpenIPSL SCRX test
+    v_0 = 1.0
+    angle_0 = 0.070492225331847
+
+    # Create GENROU machine with EFD input enabled for exciter control
+    @named genrou = PSSE_GENROU(;
+        # System parameters
+        S_b = 100e6,
+        M_b = 100e6,
+        # GENROU machine parameters from OpenIPSL SCRX.mo (lines 6-28)
+        H = 4.28,
+        D = 0,
+        Xd = 1.84,
+        Xq = 1.75,
+        Xpd = 0.41,
+        Xpq = 0.6,
+        Xppd = 0.2,
+        Xppq = 0.2,
+        Xl = 0.12,
+        Tpd0 = 5,
+        Tppd0 = 0.07,
+        Tpq0 = 0.9,
+        Tppq0 = 0.09,
+        S10 = 0.11,
+        S12 = 0.39,
+        R_a = 0,
+        pmech_input = false,
+        efd_input = true  # EFD controlled by SCRX
+    )
+
+    # Create SCRX exciter with parameters matching OpenIPSL test (lines 30-38)
+    @named scrx = PSSE_SCRX(;
+        T_AT_B = 0.1,
+        T_B = 1,
+        K = 100,
+        T_E = 0.005,
+        E_MIN = -10,
+        E_MAX = 10,
+        r_cr_fd = 0,
+        C_SWITCH = false,  # Bus fed mode
+        vothsg_input = false,  # No other signal input
+        vuel_input = false,    # No under-excitation limiter
+        voel_input = false     # No over-excitation limiter
+    )
+
+    con = [
+        # Connect exciter output to machine field voltage
+        connect(scrx.EFD_out, genrou.EFD_in)
+
+        # Connect machine field current to exciter
+        connect(genrou.XADIFD_out, scrx.XADIFD_in)
+
+        # Connect terminal voltage measurement
+        connect(genrou.ETERM_out, scrx.ECOMP_in)
+    ]
+
+    # Create bus model with proper connections
+    busmodel = MTKBus([genrou, scrx], con; name=:GEN1)
+    bm = compile_bus(busmodel, pf=pfSlack(V=v_0, δ=angle_0))
+end
+
+sol = OpenIPSL_SMIB(BUS);
+
+## Validation tests for GENROU machine variables (core 3 variables)
+@test ref_rms_error(sol, ref, VIndex(:GEN1, :genrou₊w), "gENROU.w") < 1e-5
+@test ref_rms_error(sol, ref, VIndex(:GEN1, :genrou₊delta), "gENROU.delta") < 1e-3
+@test ref_rms_error(sol, ref, VIndex(:GEN1, :genrou₊Vt), "gENROU.Vt") < 1e-4
+
+## Validation tests for SCRX exciter variables (4 control path variables)
+@test ref_rms_error(sol, ref, VIndex(:GEN1, :scrx₊leadlag₊out), "sCRX.imLeadLag.y") < 5e-5
+@test ref_rms_error(sol, ref, VIndex(:GEN1, :scrx₊amplifier₊out), "sCRX.simpleLagLim.y") < 5e-3
+@test ref_rms_error(sol, ref, VIndex(:GEN1, :genrou₊XADIFD_out₊u), "gENROU.XADIFD") < 1e-3
+@test ref_rms_error(sol, ref, VIndex(:GEN1, :scrx₊EFD_out₊u), "sCRX.EFD") < 2e-3
+
+#=
+fig_scrx = let
+    fig = Figure(size=(1400, 1200))
+    ts = refine_timeseries(sol.t)
+
+    # Plot 1: Generator Rotor Angle
+    ax1 = Axis(fig[1,1]; xlabel="Time [s]", ylabel="δ [rad]", title="Generator: Rotor Angle")
+    lines!(ax1, ref.time, ref[!, Symbol("gENROU.delta")]; label="OpenIPSL", color=Cycled(1), linewidth=2, alpha=0.7)
+    lines!(ax1, ts, sol(ts, idxs=VIndex(:GEN1, :genrou₊delta)).u; label="PowerDynamics", color=Cycled(1), linewidth=2, linestyle=:dash)
+    axislegend(ax1)
+
+    # Plot 2: Generator Speed Deviation
+    ax2 = Axis(fig[1,2]; xlabel="Time [s]", ylabel="ω [pu]", title="Generator: Speed Deviation")
+    lines!(ax2, ref.time, ref[!, Symbol("gENROU.w")]; label="OpenIPSL", color=Cycled(2), linewidth=2, alpha=0.7)
+    lines!(ax2, ts, sol(ts, idxs=VIndex(:GEN1, :genrou₊w)).u; label="PowerDynamics", color=Cycled(2), linewidth=2, linestyle=:dash)
+    axislegend(ax2)
+
+    # Plot 3: Generator Terminal Voltage
+    ax3 = Axis(fig[2,1]; xlabel="Time [s]", ylabel="Vt [pu]", title="Generator: Terminal Voltage")
+    lines!(ax3, ref.time, ref[!, Symbol("gENROU.Vt")]; label="OpenIPSL", color=Cycled(3), linewidth=2, alpha=0.7)
+    lines!(ax3, ts, sol(ts, idxs=VIndex(:GEN1, :genrou₊Vt)).u; label="PowerDynamics", color=Cycled(3), linewidth=2, linestyle=:dash)
+    axislegend(ax3)
+
+    # Plot 4: SCRX LeadLag Output
+    ax4 = Axis(fig[2,2]; xlabel="Time [s]", ylabel="[pu]", title="SCRX: LeadLag Output")
+    lines!(ax4, ref.time, ref[!, Symbol("sCRX.imLeadLag.y")]; label="OpenIPSL", color=Cycled(4), linewidth=2, alpha=0.7)
+    lines!(ax4, ts, sol(ts, idxs=VIndex(:GEN1, :scrx₊leadlag₊out)).u; label="PowerDynamics", color=Cycled(4), linewidth=2, linestyle=:dash)
+    axislegend(ax4)
+
+    # Plot 5: SCRX Amplifier Output
+    ax5 = Axis(fig[3,1]; xlabel="Time [s]", ylabel="[pu]", title="SCRX: Amplifier Output")
+    lines!(ax5, ref.time, ref[!, Symbol("sCRX.simpleLagLim.y")]; label="OpenIPSL", color=Cycled(5), linewidth=2, alpha=0.7)
+    lines!(ax5, ts, sol(ts, idxs=VIndex(:GEN1, :scrx₊amplifier₊out)).u; label="PowerDynamics", color=Cycled(5), linewidth=2, linestyle=:dash)
+    axislegend(ax5)
+
+    # Plot 6: SCRX Final EFD Output
+    ax6 = Axis(fig[3,2]; xlabel="Time [s]", ylabel="EFD [pu]", title="SCRX: Final Field Voltage")
+    lines!(ax6, ref.time, ref[!, Symbol("sCRX.EFD")]; label="OpenIPSL", color=Cycled(6), linewidth=2, alpha=0.7)
+    lines!(ax6, ts, sol(ts, idxs=VIndex(:GEN1, :scrx₊EFD_out₊u)).u; label="PowerDynamics", color=Cycled(6), linewidth=2, linestyle=:dash)
+    axislegend(ax6)
+
+    fig
+end
+=#
diff --git a/test/OpenIPSL_test/SCRX/model_simulation.mos b/test/OpenIPSL_test/SCRX/model_simulation.mos
@@ -0,0 +1,23 @@
+// Load OpenIPSL library
+loadFile(getEnvironmentVar("OPENIPSL_PATH") + "/package.mo");
+
+// Run SCRX simulation with model's native experiment settings
+simulate(OpenIPSL.Tests.Controls.PSSE.ES.SCRX);
+
+// Extract minimal machine and exciter data for testing
+filterSimulationResults(
+      "OpenIPSL.Tests.Controls.PSSE.ES.SCRX_res.mat",
+      "modelica_results.csv",
+      {
+            // GENROU machine variables (key dynamics)
+            "gENROU.delta",                 // Rotor angle
+            "gENROU.w",                     // Speed deviation (omega)
+            "gENROU.Vt",                    // Terminal voltage magnitude
+
+            // SCRX exciter variables (control path)
+            "sCRX.imLeadLag.y",            // LeadLag compensator output
+            "sCRX.simpleLagLim.y",         // Amplifier output
+            "gENROU.XADIFD",               // Field current (machine side)
+            "sCRX.EFD"                     // Final exciter field voltage output
+      }
+);
diff --git a/test/OpenIPSL_test/SCRX/modelica_results.csv.gz b/test/OpenIPSL_test/SCRX/modelica_results.csv.gz
diff --git a/test/OpenIPSL_test/SCRX/ref_data_gen.sh b/test/OpenIPSL_test/SCRX/ref_data_gen.sh
@@ -0,0 +1,38 @@
+#!/bin/bash
+
+# Get script directory
+SCRIPTDIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+
+# Set temporary directory
+TMPDIR="/tmp/scrx_sim"
+
+# Create temporary directory
+mkdir -p "$TMPDIR" && cd "$TMPDIR"
+
+# Remove existing OpenIPSL directory if it exists
+rm -rf OpenIPSL
+
+# Clone OpenIPSL at specific commit for reproducibility
+git clone --depth 1 https://github.com/OpenIPSL/OpenIPSL.git
+cd OpenIPSL
+git fetch --depth 1 origin fe8aa5c
+git checkout fe8aa5c
+cd ..
+
+# Set OpenIPSL library path
+export OPENIPSL_PATH="$TMPDIR/OpenIPSL/OpenIPSL"
+
+# Run OpenModelica simulation
+omc "$SCRIPTDIR/model_simulation.mos"
+
+# Copy results back, compress, and cleanup
+# Copy minimal data (for git repo)
+cp modelica_results.csv "$SCRIPTDIR/" || { echo "Error: Failed to copy simulation results"; exit 1; }
+
+# Remove existing compressed files if they exist to avoid conflicts
+rm -f "$SCRIPTDIR/modelica_results.csv.gz"
+
+# Compress the file
+gzip "$SCRIPTDIR/modelica_results.csv"
+
+rm -rf "$TMPDIR"
