@@ -162,7 +162,6 @@ namespace GridKit
162162 }
163163
164164 check (Ka_ > ZERO <RealT>, " Ka must be positive" );
165- check (Tc_ >= ZERO <RealT>, " Tc must be non-negative" );
166165 check (Vrmin_ <= Vrmax_, " Vrmin must be less than or equal to Vrmax" );
167166 check (UEL_ >= static_cast <IdxT>(0 ) && UEL_ <= static_cast <IdxT>(3 ),
168167 " UEL must be 0, 1, 2, or 3" );
@@ -176,14 +175,11 @@ namespace GridKit
176175 {
177176 check (E1_ > ZERO <RealT>, " E1 must be positive when saturation is enabled" );
178177 check (E2_ > ZERO <RealT>, " E2 must be positive when saturation is enabled" );
179- check (Se1_ > ZERO <RealT>, " Se1 must be positive when saturation is enabled" );
180- check (Se2_ > ZERO <RealT>, " Se2 must be positive when saturation is enabled" );
181-
182- const bool saturation_points_are_ordered =
183- (E2_ > E1_ && Se2_ > Se1_)
184- || (E2_ < E1_ && Se2_ < Se1_);
185- check (saturation_points_are_ordered,
186- " E1/E2 and Se1/Se2 must be ordered consistently" );
178+ check (Se1_ >= ZERO <RealT>, " Se1 must be non-negative when saturation is enabled" );
179+ check (Se2_ >= ZERO <RealT>, " Se2 must be non-negative when saturation is enabled" );
180+
181+ const bool sat_ordered = (E2_ > E1_ && Se2_ > Se1_) || (E2_ < E1_ && Se2_ < Se1_);
182+ check (sat_ordered, " E1/E2 and Se1/Se2 must be ordered consistently" );
187183 }
188184
189185 if (!signals_.template isAssigned <Esdc1aInternalVariables::EFD >())
@@ -314,8 +310,15 @@ namespace GridKit
314310 return 1 ;
315311 }
316312
317- const ScalarT se0 = SB_ * Math::qramp (efdp0 - SA_ );
318- const ScalarT vfe0 = (Ke_ + se0) * efdp0;
313+ const ScalarT se0 = SB_ * Math::qramp (efdp0 - SA_ );
314+
315+ RealT ke0 = Ke_;
316+ if (ke0 == ZERO <RealT>)
317+ {
318+ ke0 = (Vrmax_ / 10.0 - static_cast <RealT>(se0)) / static_cast <RealT>(efdp0);
319+ }
320+
321+ const ScalarT vfe0 = ke0 * efdp0 + se0;
319322 const ScalarT vr0 = vfe0;
320323 const ScalarT vhv0 = vr0 / Ka_;
321324
@@ -344,15 +347,16 @@ namespace GridKit
344347 vll0 = vuel0 + inverseRamp (gate_margin0);
345348 }
346349
347- const ScalarT vf0 = ScalarT{ZERO <RealT>};
348350 const ScalarT ev0 = vll0;
349351 const ScalarT xll0 = ev0;
350- const ScalarT vref0 = ev0 + vc0 + vf0 - vs0 - uel_on_ * vuel0;
352+ const ScalarT vref0 = ev0 + vc0 - vs0 - uel_on_ * vuel0;
353+
354+ Ke_ = ke0;
351355
352356 y[EFDP ] = efdp0;
353357 y[VC ] = vc0;
354358 y[VR ] = vr0;
355- y[VF ] = vf0 ;
359+ y[VF ] = ZERO <RealT> ;
356360 y[XLL ] = xll0;
357361 y[EV ] = ev0;
358362 y[VLL ] = vll0;
@@ -560,7 +564,7 @@ namespace GridKit
560564
561565 const ScalarT ec = std::sqrt (wb[0 ] * wb[0 ] + wb[1 ] * wb[1 ]);
562566 const ScalarT ev_target = vref + vs + uel_on_ * vuel - vc - vf;
563- const ScalarT vfe_target = ( Ke_ + se) * efdp;
567+ const ScalarT vfe_target = Ke_ * efdp + se ;
564568 const ScalarT efdp_rate = (vr - vfe) / Te_;
565569 const ScalarT limited_efdp_rate = awmin (efdp, efdp_rate, ZERO <RealT>);
566570
@@ -749,7 +753,7 @@ namespace GridKit
749753 * @brief Resolve the parameter-derived constants and selector masks
750754 *
751755 * Raises the transducer, regulator, lead-lag, exciter, and feedback
752- * lags to the well-posedness floor, fits the quadratic saturation
756+ * lags to the well-posedness floor, fits the scaled- quadratic saturation
753757 * curve, and turns the three selectors into multiplicative masks. The
754758 * masks let the residual select signal routing without
755759 * parameter-dependent control flow, which keeps its structure fixed for
@@ -814,23 +818,37 @@ namespace GridKit
814818 // The disabled test matches the verify() predicate exactly.
815819 const bool saturation_disabled =
816820 Se1_ == ZERO <RealT> && Se2_ == ZERO <RealT>;
817- const bool saturation_points_are_ordered =
818- (E2_ > E1_ && Se2_ > Se1_)
819- || (E2_ < E1_ && Se2_ < Se1_);
821+ const bool sat_ordered = (E2_ > E1_ && Se2_ > Se1_) || (E2_ < E1_ && Se2_ < Se1_);
820822 const bool saturation_consistent =
821823 E1_ > ZERO <RealT> && E2_ > ZERO <RealT>
822- && Se1_ > ZERO <RealT> && Se2_ > ZERO <RealT>
823- && saturation_points_are_ordered ;
824+ && Se1_ >= ZERO <RealT> && Se2_ >= ZERO <RealT>
825+ && sat_ordered ;
824826 if (saturation_disabled || !saturation_consistent)
825827 {
826828 SA_ = ZERO <RealT>;
827829 SB_ = ZERO <RealT>;
828830 return ;
829831 }
830832
831- const RealT C = std::sqrt (Se2_ / Se1_);
833+ if (Se1_ == ZERO <RealT>)
834+ {
835+ const RealT dE = E2_ - E1_ ;
836+ SA_ = E1_ ;
837+ SB_ = Se2_ * E2_ / (dE * dE);
838+ return ;
839+ }
840+
841+ if (Se2_ == ZERO <RealT>)
842+ {
843+ const RealT dE = E1_ - E2_ ;
844+ SA_ = E2_ ;
845+ SB_ = Se1_ * E1_ / (dE * dE);
846+ return ;
847+ }
848+
849+ const RealT C = std::sqrt (Se2_ * E2_ / (Se1_ * E1_ ));
832850 SA_ = (C * E1_ - E2_ ) / (C - ONE <RealT>);
833- SB_ = Se1_ / ((E1_ - SA_ ) * (E1_ - SA_ ));
851+ SB_ = Se1_ * E1_ / ((E1_ - SA_ ) * (E1_ - SA_ ));
834852 }
835853
836854 /* *
@@ -840,30 +858,17 @@ namespace GridKit
840858 * *input*, so the residual reproduces the requested output through the
841859 * same smooth ramp it evaluates.
842860 *
843- * For large positive values, the ramp is effectively equal to the input, so the
844- * inverse is effectively the output. In that regime this function returns `ramp_output` directly.
845- * This branching is numerically more robust.
846- *
847861 * @param[in] ramp_output Strictly positive requested ramp output.
848862 * @return The input the smooth ramp maps to the requested output.
849863 *
850864 * @pre @p ramp_output is finite and strictly positive.
851- * @warning This function contains conditional branching and may be used
852- * during initialization, but not during residual or Jacobian
853- * evaluation.
854865 */
855866 template <typename scalar_type, typename index_type>
856867 typename Esdc1a<scalar_type, index_type>::RealT
857868 Esdc1a<scalar_type, index_type>::inverseRamp(RealT ramp_output) const
858869 {
859- static constexpr RealT SOFTPLUS_WIDTH = static_cast <RealT>(50.0 );
860-
861- const RealT scaled_output = Math::MU <RealT> * ramp_output;
862- if (scaled_output > SOFTPLUS_WIDTH )
863- {
864- return ramp_output;
865- }
866- return std::log (std::expm1 (scaled_output)) / Math::MU <RealT>;
870+ const RealT mu = Math::MU <RealT>;
871+ return ramp_output + std::log (-std::expm1 (-mu * ramp_output)) / mu;
867872 }
868873
869874 /* *
0 commit comments