1313#include < GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp>
1414#include < GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp>
1515#include < GridKit/Model/VariableMonitorImpl.hpp>
16+ #include < GridKit/Testing/Testing.hpp>
1617#include < GridKit/Utilities/Logger/Logger.hpp>
1718
1819namespace GridKit
@@ -22,6 +23,7 @@ namespace GridKit
2223 namespace Stabilizer
2324 {
2425 using Log = ::GridKit::Utilities::Logger;
26+ using Testing::isEqual;
2527
2628 template <typename scalar_type, typename index_type>
2729 Ieeest<scalar_type, index_type>::Ieeest()
@@ -191,7 +193,7 @@ namespace GridKit
191193 ret += 1 ;
192194 }
193195
194- if (a2_ == ZERO <RealT> && a3_ == ZERO <RealT> && a4_ == ZERO <RealT> && a1_ != ZERO <RealT>)
196+ if (isEqual ( a2_, ZERO <RealT>) && isEqual ( a3_, ZERO <RealT>) && isEqual ( a4_, ZERO <RealT>) && ! isEqual ( a1_, ZERO <RealT>) )
195197 {
196198 Log::error () << " Ieeest: a2, a3, and a4 are all zero - no valid notch filter\n " ;
197199 ret += 1 ;
@@ -228,7 +230,7 @@ namespace GridKit
228230 const ScalarT x6 = v5;
229231 const ScalarT v6 = v5;
230232 const ScalarT x7 = v6;
231- const ScalarT v7 = (T6_ == ZERO <RealT>) ? Ks_ * v6 : zero;
233+ const ScalarT v7 = isEqual (T6_ , ZERO <RealT>) ? Ks_ * v6 : zero;
232234
233235 y_[0 ] = x1;
234236 y_[1 ] = x2;
@@ -249,13 +251,13 @@ namespace GridKit
249251 template <typename scalar_type, typename index_type>
250252 int Ieeest<scalar_type, index_type>::tagDifferentiable()
251253 {
252- tag_[0 ] = (a2_ != ZERO <RealT> || a3_ != ZERO <RealT> || a4_ != ZERO <RealT>);
254+ tag_[0 ] = (! isEqual ( a2_, ZERO <RealT>) || ! isEqual ( a3_, ZERO <RealT>) || ! isEqual ( a4_, ZERO <RealT>) );
253255 tag_[1 ] = tag_[0 ];
254- tag_[2 ] = (a3_ != ZERO <RealT> || a4_ != ZERO <RealT>);
255- tag_[3 ] = (a4_ != ZERO <RealT>);
256- tag_[4 ] = (T2_ != ZERO <RealT>);
257- tag_[5 ] = (T4_ != ZERO <RealT>);
258- tag_[6 ] = (T6_ != ZERO <RealT>);
256+ tag_[2 ] = (! isEqual ( a3_, ZERO <RealT>) || ! isEqual ( a4_, ZERO <RealT>) );
257+ tag_[3 ] = (! isEqual ( a4_, ZERO <RealT>) );
258+ tag_[4 ] = (! isEqual ( T2_ , ZERO <RealT>) );
259+ tag_[5 ] = (! isEqual ( T4_ , ZERO <RealT>) );
260+ tag_[6 ] = (! isEqual ( T6_ , ZERO <RealT>) );
259261 tag_[7 ] = false ;
260262 tag_[8 ] = false ;
261263 tag_[9 ] = false ;
@@ -315,17 +317,26 @@ namespace GridKit
315317
316318 ScalarT u = ws[0 ];
317319
318- f[0 ] = -tag_[0 ] * x1_dot + x2;
319- f[1 ] = -tag_[1 ] * x2_dot + x3;
320- f[2 ] = -tag_[2 ] * x3_dot + x4;
320+ const ScalarT zero{ZERO <RealT>};
321+ const ScalarT one{ONE <RealT>};
322+ const ScalarT dx1 = tag_[0 ] ? one : zero;
323+ const ScalarT dx2 = tag_[1 ] ? one : zero;
324+ const ScalarT dx3 = tag_[2 ] ? one : zero;
325+ const ScalarT t2 = tag_[4 ] ? one : zero;
326+ const ScalarT t4 = tag_[5 ] ? one : zero;
327+ const ScalarT t6 = tag_[6 ] ? one : zero;
328+
329+ f[0 ] = -dx1 * x1_dot + x2;
330+ f[1 ] = -dx2 * x2_dot + x3;
331+ f[2 ] = -dx3 * x3_dot + x4;
321332 f[3 ] = -a4_ * x4_dot - x1 - a1_ * x2 - a2_ * x3 - a3_ * x4 + u;
322333 f[4 ] = -T2_ * x5_dot - x5 + v4;
323334 f[5 ] = -T4_ * x6_dot - x6 + v5;
324335 f[6 ] = -T6_ * x7_dot - x7 + v6;
325336 f[7 ] = -v4 + x1 + A5_ * x2 + A6_ * x3;
326- f[8 ] = tag_[ 4 ] * (-T2_ * (v5 - x5) + T1_ * (v4 - x5)) + (1 - tag_[ 4 ] ) * (v4 - v5);
327- f[9 ] = tag_[ 5 ] * (-T4_ * (v6 - x6) + T3_ * (v5 - x6)) + (1 - tag_[ 5 ] ) * (v5 - v6);
328- f[10 ] = tag_[ 6 ] * (-T6_ * v7 + Ks_ * T5_ * (v6 - x7)) + (1 - tag_[ 6 ] ) * (Ks_ * v6 - v7);
337+ f[8 ] = t2 * (-T2_ * (v5 - x5) + T1_ * (v4 - x5)) + (one - t2 ) * (v4 - v5);
338+ f[9 ] = t4 * (-T4_ * (v6 - x6) + T3_ * (v5 - x6)) + (one - t4 ) * (v5 - v6);
339+ f[10 ] = t6 * (-T6_ * v7 + Ks_ * T5_ * (v6 - x7)) + (one - t6 ) * (Ks_ * v6 - v7);
329340 f[11 ] = -vss + Math::clamp (v7, Lsmin_, Lsmax_);
330341
331342 return 0 ;
0 commit comments