BUG: Renewable Controller edge case

GridKit/Model/PhasorDynamics/Controller/REECB/README.md

@@ -152,7 +152,6 @@ $e_Q$                        | [p.u.]   | Reactive-power error
 $V_Q^\mathrm{PI}$           | [p.u.]   | Reactive-power PI output                      |
 $e_V^\mathrm{PI}$           | [p.u.]   | Voltage-control error                         |
 $r_P^\mathrm{ord}$          | [p.u./s] | Limited active-power order rate               | Component base
-$I_L^{\max}$                | [p.u.]   | Current-circle continuation state             | Component base
 $I_L^\mathrm{cap}$          | [p.u.]   | Off-axis current capacity                     | Component base
 $I_q^{\max}$                | [p.u.]   | Reactive-current limit                        | Component base
 $I_p^{\max}$                | [p.u.]   | Active-current limit                          | Component base
@@ -209,8 +208,7 @@ $P^\mathrm{ref}$       | [p.u.] | Unknown | External active-power reference
   0 &= -V_Q^\mathrm{PI}+\text{clamp}(K_\mathrm{qp}e_Q+x_Q^\mathrm{PI};\,V^{\min},V^{\max}) \\
   0 &= -e_V^\mathrm{PI}+s_Q^\mathrm{PI}V_Q^\mathrm{PI}+s_V^\mathrm{ref}Q^\mathrm{ext}-s_QV^\mathrm{meas} \\
   0 &= -r_P^\mathrm{ord}+\text{aslew}\!\left(\dfrac{k_\mathrm{base}P^\mathrm{ref}-P^\mathrm{ord}}{T_\mathrm{pord}};\,R_P^{\min},R_P^{\max}\right) \\
-  0 &= -I_L^{\max}\sqrt{(I_L^{\max})^2+\epsilon_0}+(I^{\max}-I^\mathrm{high})(I^{\max}+I^\mathrm{high}) \\
-  0 &= -I_L^\mathrm{cap}+\dfrac{(I_L^{\max})^2}{\sqrt{(I_L^{\max})^2+\epsilon_0}} \\
+  0 &= -I_L^\mathrm{cap}+\text{sqrtramp}\big((I^{\max}-I^\mathrm{high})(I^{\max}+I^\mathrm{high})\big) \\
   0 &= -I_q^{\max}+s_\mathrm{pq}I_L^\mathrm{cap}+s_\mathrm{pq}^\mathrm{off}I^{\max} \\
   0 &= -I_p^{\max}+s_\mathrm{pq}I^{\max}+s_\mathrm{pq}^\mathrm{off}I_L^\mathrm{cap} \\
   0 &= -I_q^\mathrm{base}+\text{clamp}(K_\mathrm{vp}e_V^\mathrm{PI}+x_V^\mathrm{PI};\,-I_q^{\max},I_q^{\max}) \\
@@ -224,6 +222,7 @@ Here
 $I^\mathrm{high}=s_\mathrm{pq}k_\mathrm{base}I_p^\mathrm{cmd}
 +s_\mathrm{pq}^\mathrm{off}k_\mathrm{base}I_q^\mathrm{cmd}$ and
 $\epsilon_0=100\epsilon_\mathrm{machine}$.
+[Appendix C](#appendix-c-sqrtramp) defines `sqrtramp`.
 
 CommonMath defines the [`antiwindup`](../../../../CommonMath.md#antiwindup) and
 [smooth limiter](../../../../CommonMath.md#derived-functions) functions used in
@@ -272,23 +271,23 @@ and $I_q=k_\mathrm{base}I_q^\mathrm{cmd}$ be the component-base commands.
 
 Let $I_q^\mathrm{need}$ be the smooth-clamp input magnitude required to
 reproduce $I_q$, including the Volt/VAr recovery margin. Initialization raises
-$I^\max$, if needed, to the smallest finite limit satisfying
+$I^\max$, if needed, to a finite limit satisfying
 
 ```math
 \begin{aligned}
   h &= s_\mathrm{pq}I_p+s_\mathrm{pq}^\mathrm{off}|I_q| \\
   \ell &= s_\mathrm{pq}I_q^\mathrm{need}+s_\mathrm{pq}^\mathrm{off}I_p \\
   I &\ge \text{max}(I^\max,h,\ell,I_q^\mathrm{need}) \\
-  x\sqrt{x^2+\epsilon_0} &= (I-h)(I+h) \\
-  \dfrac{x^2}{\sqrt{x^2+\epsilon_0}} &\ge \ell.
+  \text{sqrtramp}\big((I-h)(I+h)\big) &\ge \ell.
 \end{aligned}
 ```
 
-The solution initializes $I^\max\leftarrow I$, $I_L^\max\leftarrow x$, and
-$I_L^\mathrm{cap}\leftarrow x^2/\sqrt{x^2+\epsilon_0}$. The operating point is
-rejected if no finite solution exists or the reactive-current injection is
-incompatible. Q, V, and P limits are expanded as needed; each adjustment logs
-a warning.
+The solution initializes
+$I^\max\leftarrow I$ and
+$I_L^\mathrm{cap}\leftarrow\text{sqrtramp}\big((I-h)(I+h)\big)$. The
+operating point is rejected if no finite solution exists or the
+reactive-current injection is incompatible. Q, V, and P limits are expanded as
+needed; each adjustment logs a warning.
 
 ```math
 \begin{aligned}
@@ -375,11 +374,11 @@ Output  | Units  | Description                     | Note
 - `initializationAndSignals()` checks initialization, signals, monitors, and power bases.
 - `initializationDomain()` checks rejected inputs and limit expansion.
 - `initializationExactness()` checks endpoint and current-circle initialization.
-- `residualEquations()` checks the fixed residual answer key.
+- `residualEquations()` checks residual equations at a hand-computable state.
 - `selectorConfigurations()` checks selectors and optional ports.
 - `voltVarReferenceBase()` checks `qext` units.
 - `reactiveControl()` checks the reactive-control paths.
-- `activeCurrentControl()` checks active-current control and current priority.
+- `activeCurrentControl()` checks active-current control, current priority, and the current circle.
 - `jacobian()` compares the Enzyme and dependency-tracking Jacobians.
 - `regcaReecbRepca()` checks the public-signal REGCA-REECB-REPCA control loop.
 - `reecb()` checks construction through the production system-data path.
@@ -412,6 +411,26 @@ where $\rho$ is GridKit's smooth
 reduces to $\text{clamp}(f;\ell,u)$; the smooth form preserves
 $\text{aslew}(0;\ell,u)=0$.
 
+## Appendix C: `sqrtramp`
+
+The off-axis capacity approaches $\sqrt{\max(0,s)}$, where
+$s=(I^{\max}-I^\mathrm{high})(I^{\max}+I^\mathrm{high})$. Let
+$\kappa=\epsilon_0/\mu$, with $\mu$ GridKit's
+[smoothing scale](../../../../CommonMath.md#primitives). REECB uses
+
+```math
+\begin{aligned}
+  h_\kappa(s) &= \tfrac{1}{2}\left(s+\sqrt{s^2+\kappa^2}\right), \\
+  \text{sqrtramp}(s) &=
+    \dfrac{\mu h_\kappa(s)}{\sqrt{\mu^2 h_\kappa(s)+1}}.
+\end{aligned}
+```
+
+Both softenings are differentiable for every finite $s$. The narrow hinge
+leaves a negligible positive capacity when the circle is closed or
+over-driven, while preserving the current-command dependence through the
+transition. Increasing $\mu$ recovers the exact current-circle leg.
+
 [^wecc-reecb-specification]: [WECC REMTF, *Generic Solar Photovoltaic System Dynamic Simulation Model Specification*](https://www.wecc.org/sites/default/files/documents/meeting/2024/WECC-Solar-PV-Dynamic-Model-Specification-September-2012.pdf), September 2012.
 
 [^reecb-mva-base]: The [WECC Central Station Photovoltaic Power Plant Model Validation Guideline](https://www.wecc.org/sites/default/files/documents/program/2024/Central%20Station%20Photovoltaic%20Power%20Plant%20Model%20Validation%20Guideline%20June%2017%202015.pdf)

GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp

@@ -44,7 +44,6 @@ namespace GridKit
         VPIQ,   ///< \f$V_Q^\mathrm{PI}\f$ Algebraic reactive-power PI output [p.u.]
         EPIV,   ///< \f$e_V^\mathrm{PI}\f$ Algebraic voltage-control error [p.u.]
         RPORD,  ///< \f$r_P^\mathrm{ord}\f$ Algebraic limited active-power order rate [p.u./s]
-        ILMAX,  ///< \f$I_L^\max\f$ Algebraic current-circle continuation state on component base [p.u.]
         ILCAP,  ///< \f$I_L^\mathrm{cap}\f$ Algebraic off-axis current capacity on component base [p.u.]
         IQMAX,  ///< \f$I_q^\max\f$ Algebraic reactive-current limit on component base [p.u.]
         IPMAX,  ///< \f$I_p^\max\f$ Algebraic active-current limit on component base [p.u.]
@@ -106,7 +105,7 @@ namespace GridKit
         using InternalVariablesT = ReecbInternalVariables;
         using ExternalVariablesT = ReecbExternalVariables;
 
-        /// Current-circle regularization and initialization reconstruction tolerance.
+        /// Initialization reconstruction tolerance.
         static constexpr RealT INITIALIZATION_TOLERANCE =
             static_cast<RealT>(100.0) * std::numeric_limits<RealT>::epsilon();
 
@@ -139,12 +138,16 @@ namespace GridKit
             ScalarT*       f);
 
       private:
+        /// Hinge width chosen so closed-circle leakage stays below the
+        /// initialization tolerance [p.u. current squared].
+        static constexpr RealT CURRENT_CIRCLE_KNEE =
+            INITIALIZATION_TOLERANCE / Math::MU<RealT>;
+
         struct InitialPoint;
 
         struct InitialCurrentLimit
         {
           RealT total_limit;
-          RealT continuation;
           RealT off_axis_capacity;
         };
 
@@ -154,17 +157,18 @@ namespace GridKit
         /// Smooth anti-windup derivative within a moving symmetric band.
         [[gnu::always_inline]] static inline ScalarT awband(ScalarT state, ScalarT rate, ScalarT band);
 
-        /// Current-circle continuation state for an initial component-base limit.
-        static RealT circleState(RealT imax, RealT high);
+        /// Smooth nonnegative root of a current-circle squared radius.
+        template <typename ValueT>
+        [[gnu::always_inline]] static inline ValueT sqrtramp(ValueT x);
 
-        /// Off-axis component-base capacity provided by a continuation state.
-        static RealT capacity(RealT ilmax);
+        /// Analytic nonnegative seed for inverting sqrtramp().
+        static RealT isqrtramp(RealT y);
 
-        /// Bisect an initial-limit bracket to its first upper-side point.
-        template <typename FuncT>
-        static RealT bisect(RealT a, RealT b, FuncT below);
+        /// Overflow-resistant difference of squares.
+        template <typename ValueT>
+        [[gnu::always_inline]] static inline ValueT circleSquare(RealT limit, ValueT high);
 
-        /// Solve the smallest feasible initial limit at or above `lower`.
+        /// Solve a feasible initial limit at or above `lower`.
         static std::optional<InitialCurrentLimit> solveInitialLimit(RealT lower, RealT high, RealT low);
 
         bool buildInitialPoint(InitialPoint& point);