Add split-SSN GFL model and state-space extraction for DP Ph3 - #698
Add split-SSN GFL model and state-space extraction for DP Ph3#698georgii-tishenin wants to merge 2 commits into
Conversation
Signed-off-by: Georgii Tishenin <georgii.tishenin@eonerc.rwth-aachen.de>
Signed-off-by: Georgii Tishenin <georgii.tishenin@eonerc.rwth-aachen.de>
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
There was a problem hiding this comment.
DPsim LLM review
Claim vs. code: The diff appears to implement the claimed feature set, but the provided patch is truncated so the example’s full validation/comparison behavior cannot be fully verified from the diff shown.
TL;DR: Three real issues stand out: unchecked numerical solves and parameter divisions in the split-SSN base, plus per-step heap allocations in the hot simulation path; the remaining items are lower-confidence documentation/binding gaps or duplicate variants of the same concerns. The state-space and stamping changes otherwise surfaced no separate correctness findings beyond these numerics and real-time determinism risks.
Found 6 high, 4 medium, 1 low (8 anchored to lines below).
🔴 Critical & high
- Per-step dynamic allocations in MNA post‑step
[high · 80% confidence]indpsim-models/src/DP/DP_VTypeSplitSSNComp.cpp:294(details inline) - Per-step dynamic allocations in MNA pre‑step
[high · 77% confidence]indpsim-models/src/DP/DP_VTypeSplitSSNComp.cpp:267(details inline) - Avoid dynamic memory allocation in the hot path of mnaCompPreStep
[high · 76% confidence]indpsim-models/src/DP/DP_VTypeSplitSSNComp.cpp:267(details inline) - Avoid dynamic memory allocation in the hot path of mnaCompPostStep
[high · 74% confidence]indpsim-models/src/DP/DP_VTypeSplitSSNComp.cpp:299(details inline) - Division by mRc without finiteness guard
[high · 72% confidence]indpsim-models/src/DP/DP_Ph3_SSN_GFL_Split.cpp:329(details inline) - Use DOUBLE_EPSILON instead of hard-coded 1e-9 tolerance
[high · 68% confidence]indpsim-models/src/DP/DP_Ph3_SSN_GFL_Split.cpp:303(details inline)
🟡 Suggestions
- Check LU factorization success and handle rank-deficient cases
[medium · 82% confidence]indpsim-models/src/DP/DP_VTypeSplitSSNComp.cpp:205(details inline) - Potential division by zero in mdBControllerOutput calculation
[medium · 80% confidence]indpsim-models/src/DP/DP_VTypeSplitSSNComp.cpp:205(details inline) - Potential division by zero in mdBControllerOutput solve
[medium · 78% confidence]indpsim-models/src/DP/DP_VTypeSplitSSNComp.cpp:205
The line mdBControllerOutput = lhs.fullPivLu().solve(0.5 * mTimeStep * mBControllerOutput) can fail or produce NaN if lhs is singular or ill-conditioned, especially if timeStep is zero or parameters yield a singular plant matrix. No finiteness guard is present.
Fix: Add a check that mTimeStep > DOUBLE_EPSILON and that lhs is well-conditioned (e.g., its smallest singular value exceeds DOUBLE_EPSILON). If not, throw std::invalid_argument with a descriptive message. - Potential numerical instability in state update due to unchecked matrix inversion
[medium · 76% confidence]indpsim-models/src/DP/DP_VTypeSplitSSNComp.cpp:205
The line**mX = mdA * (**mX) + mdB * (uNew + uOld) + 2.0 * mdBControllerOutput * delayedOutputUsed;updates the state using a matrix inversion (via mdBControllerOutput) without checking the condition number or invertibility of the matrixlhs = identity - 0.5 * mTimeStep * mAinrecomputeDiscreteModel(). If the matrix is ill-conditioned or singular, the state update can become unstable or incorrect.
Fix: Add a check afterlhs.fullPivLu().solve(...)to ensure the matrix is invertible (e.g., check the condition number or the success of the solve). If the matrix is singular, throw an exception or use a fallback (e.g., pseudo-inverse). - No singularity check for lhs matrix solve
[low · 79% confidence]indpsim-models/src/DP/DP_VTypeSplitSSNComp.cpp:205
recomputeDiscreteModel() solves lhs.fullPivLu().solve(...). If lhs is singular (e.g., due to extreme parameter values or zero time step), the solution may contain NaNs without detection, silently corrupting the discrete model.
Fix: After solving, verify that the solution is finite, e.g., using Math::isFinite on mdBControllerOutput, and raise an exception if not.
Not shown: 6 tentative. All are in the workflow log.
Claim vs. implementation
- Claimed: Add split-SSN support for three-phase DP voltage-type components and a DP Ph3 split grid-following inverter, plus state-space extraction and an example.
- Done: Adds a new split-SSN base and three-phase two-terminal specialization, implements DP Ph3 SSN_GFL_Split with delayed controller output and analytical controller/state-space hooks, and adds a DP Ph3 state-space extraction example.
- Difference: The diff appears to implement the claimed feature set, but the provided patch is truncated so the example’s full validation/comparison behavior cannot be fully verified from the diff shown.
How this review was produced
13 specialized finder passes raised 107 findings over the diff and the full changed sources. After de-duplication, 107 were re-checked against the current file and the base-class / interface headers it inherits (code as truth), escalating survivors to a stronger model: 90 refuted as unsupported, 17 kept (6 tentative).
Refuted by verification:
- Incorrect sign in PLL state derivative (dpsim-models/src/DP/DP_Ph3_SSN_GFL_Split.cpp): f(PhiPLL, 0) is explicitly set to vcDq.imag() in the source, so the claim of an inverted sign is contradicted by the file
- Incorrect power-filter dynamics (dpsim-models/src/DP/DP_Ph3_SSN_GFL_Split.cpp): the code implements f(PFiltered, 0) = mOmegaCutoff * (pInst - x(PFiltered, 0)), which is the standard first-order low-pass form
- Incorrect current measurement projection in controller output (dpsim-models/src/DP/DP_Ph3_SSN_GFL_Split.cpp): evaluateControllerOutput uses the same rot/projection pattern as the derivative path and the code already conjugates projection[p] when mapping back to phase outputs
- Incorrect rotation sign in analytical Jacobian (dpsim-models/src/DP/DP_Ph3_SSN_GFL_Split.cpp): calculateControllerAnalyticalJacobians uses std::exp(-j * x(Psi, 0)), which is the same rotation convention as the other controller methods
- Incorrect PLL initial frequency error (dpsim-models/src/DP/DP_Ph3_SSN_GFL_Split.cpp): the initialization of PhiPLL to (omega - mOmegaN) / mKiPLL is present in the file and matches the controller’s PLL state convention used here
- ... and 10 more, in the workflow log.
Automated, non-blocking review. May be wrong. Models: find mistral-small-4-119b-2603, gpt-oss-120b → verify gpt-5.4-mini → final gpt-5.5.
| mnaCompUpdateVoltage(**leftVector); | ||
| mnaCompUpdateCurrent(**leftVector); | ||
|
|
||
| const Matrix uOld = packComplex(uOldComplex); |
There was a problem hiding this comment.
Per-step dynamic allocations in MNA post‑step
severity: high · confidence: 80%
The method DP::VTypeSplitSSNComp::mnaCompPostStep allocates multiple temporary matrices each step (uOld, uNew, measurementNew, derivative, output, etc.) and calls functions that return newly allocated matrices (buildControllerMeasurement, evaluateControllerStateDerivative, evaluateControllerOutput). This introduces heap allocations in the critical simulation loop, harming real‑time determinism.
Suggested fix: Introduce reusable member matrices for uOld, uNew, measurementNew, derivative, and output. Refactor buildControllerMeasurement, evaluateControllerStateDerivative, and evaluateControllerOutput to write into pre‑allocated output arguments instead of returning new matrices.
Checked against the source: mnaCompPostStep creates uOld, uNew, measurementNew and other dynamic Matrix temporaries each step
stage: realtime-resources
| } | ||
|
|
||
| void DP::VTypeSplitSSNComp::mnaCompPreStep(Real, Int) { | ||
| const Matrix history = |
There was a problem hiding this comment.
Per-step dynamic allocations in MNA pre‑step
severity: high · confidence: 77%
The method DP::VTypeSplitSSNComp::mnaCompPreStep creates several temporary Eigen matrices (history, result of packComplex, etc.) on each simulation step. These allocations occur inside the hot path of the DP solver and can cause non‑deterministic timing, breaking real‑time simulation guarantees.
Suggested fix: Pre‑allocate member buffers for the intermediate results (e.g., a Matrix for history, a Matrix for the packed voltage, and a Matrix for the controller measurement) and reuse them each step. Prefer fixed‑size Eigen types (e.g., Eigen::Matrix<Real,12,1>) when dimensions are known at compile time to avoid heap allocations.
Checked against the source: the pre-step hot path allocates Matrix history and packComplex/matrix-product results
stage: realtime-resources
| } | ||
|
|
||
| void DP::VTypeSplitSSNComp::mnaCompPreStep(Real, Int) { | ||
| const Matrix history = |
There was a problem hiding this comment.
Avoid dynamic memory allocation in the hot path of mnaCompPreStep
severity: high · confidence: 76%
In VTypeSplitSSNComp::mnaCompPreStep, the expression mC * (mdA * (**mX) + mdB * packComplex(**mIntfVoltage) + 2.0 * mdBControllerOutput * mDelayedControllerOutput) builds multiple temporaries and may allocate memory for the result of packComplex and the matrix-vector products. This breaks real-time determinism.
Suggested fix: Preallocate a temporary for the plant+delay state term and compute it in-place. Avoid packComplex inside the step loop by working with packed-real views of **mIntfVoltage and **mIntfCurrent. Consider fusing the controller-output history term into a single GEMV call.
Checked against the source: mnaCompPreStep creates a Matrix history from nested matrix products and packComplex() every step
stage: realtime-resources
| const Matrix delayedOutputUsed = mDelayedControllerOutput; | ||
| const Matrix controllerBdUsed = mControllerBd; | ||
|
|
||
| **mX = mdA * (**mX) + mdB * (uNew + uOld) + |
There was a problem hiding this comment.
Avoid dynamic memory allocation in the hot path of mnaCompPostStep
severity: high · confidence: 74%
In VTypeSplitSSNComp::mnaCompPostStep, the expression mdA * (**mX) + mdB * (uNew + uOld) + 2.0 * mdBControllerOutput * delayedOutputUsed performs a temporary copy of (**mX) and uNew+uOld, which can trigger dynamic memory allocation for large state vectors. This breaks real-time determinism.
Suggested fix: Rewrite the state update to compute in-place into a preallocated temporary or directly into **mX, avoiding temporaries. Use noalias() or expression templates if Eigen is configured accordingly, or compute the update in a single AXPY-style call.
Checked against the source: mnaCompPostStep constructs dynamic Eigen temporaries in the state-update expression each step
stage: realtime-resources
| (**mX)(mIfIm[p], 0) = ifAbc(p, 0).imag(); | ||
| } | ||
| **mIntfVoltage = u; | ||
| **mIntfCurrent = (u - vcAbc) / mRc; |
There was a problem hiding this comment.
Division by mRc without finiteness guard
severity: high · confidence: 72%
The interface current is computed as (u - vcAbc) / mRc. If mRc is zero or near-zero due to invalid parameters, this division can produce NaN or Inf, corrupting the admittance matrix and violating numerical soundness. DPsim guidelines require defensive validation and finiteness checks for any division involving model parameters.
Suggested fix: Add a runtime check in setParameters to ensure mRc > DOUBLE_EPSILON. If violated, throw std::invalid_argument with a clear message. Also add Math::isFinite checks before the division in initializeFromNodesAndTerminals.
Checked against the source: mIntfCurrent divides by mRc and setParameters only checks rc <= 0.0, not near-zero or finite values
stage: numerics
| const MatrixComp u = buildInitialInputFromNodes(frequency); | ||
| Complex vc = u(0, 0), iGrid(0.0, 0.0); | ||
| for (Int iteration = 0; iteration < 10; ++iteration) { | ||
| if (std::abs(vc) < 1e-9) |
There was a problem hiding this comment.
Use DOUBLE_EPSILON instead of hard-coded 1e-9 tolerance
severity: high · confidence: 68%
The initialization loop in initializeFromNodesAndTerminals uses a hard-coded tolerance of 1e-9 to check convergence of the operating-point iteration. DPsim conventions require using the shared constant DOUBLE_EPSILON from Definitions.h for all near-zero and tolerance comparisons to ensure consistency and avoid magic numbers.
Suggested fix: Replace the literal 1e-9 with DOUBLE_EPSILON in the convergence check at line 308. Also replace the similar literal 1e-9 at line 303 with DOUBLE_EPSILON for symmetry.
Checked against the source: initializeFromNodesAndTerminals uses the hard-coded tolerance 1e-9 instead of the shared DOUBLE_EPSILON constant
stage: numerics
| const Matrix identity = Matrix::Identity(mA.rows(), mA.cols()); | ||
| const Matrix lhs = identity - 0.5 * mTimeStep * mA; | ||
| mdBControllerOutput = | ||
| lhs.fullPivLu().solve(0.5 * mTimeStep * mBControllerOutput); |
There was a problem hiding this comment.
Check LU factorization success and handle rank-deficient cases
severity: medium · confidence: 82%
The call to fullPivLu().solve(...) does not check the success of the factorization. If the plant matrix A is ill-conditioned or singular, the solver may return NaNs or silently incorrect values, corrupting the simulation state without warning.
Suggested fix: Check solver.info() after solve and throw a descriptive exception on failure. Consider falling back to a more robust solver (e.g., pseudo-inverse) if the matrix is singular, but always raise an error in debug builds.
Checked against the source: the FullPivLU solve result is used directly with no success/rank or NaN check
stage: realtime-resources
| const Matrix identity = Matrix::Identity(mA.rows(), mA.cols()); | ||
| const Matrix lhs = identity - 0.5 * mTimeStep * mA; | ||
| mdBControllerOutput = | ||
| lhs.fullPivLu().solve(0.5 * mTimeStep * mBControllerOutput); |
There was a problem hiding this comment.
Potential division by zero in mdBControllerOutput calculation
severity: medium · confidence: 80%
The line mdBControllerOutput = lhs.fullPivLu().solve(0.5 * mTimeStep * mBControllerOutput); does not check if lhs is invertible or if the solve succeeds. If lhs is singular or nearly singular, the result can be NaN or incorrect, leading to instability.
Suggested fix: Check the condition number of lhs or the success of the solve. If the matrix is singular, throw an exception or use a fallback (e.g., pseudo-inverse).
Checked against the source: lhs.fullPivLu().solve(...) has no invertibility or result-finiteness guard
stage: io-robustness
|
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #698 +/- ##
==========================================
- Coverage 72.59% 71.43% -1.17%
==========================================
Files 522 526 +4
Lines 36118 36722 +604
Branches 19785 19978 +193
==========================================
+ Hits 26221 26232 +11
- Misses 9812 10414 +602
+ Partials 85 76 -9 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
1 participant
This PR adds split-SSN support for three-phase dynamic-phasor voltage-type components and introduces a corresponding DP Ph3 grid-following inverter model.
The new split formulation separates the fixed linear electrical plant from the nonlinear controller. The controller is linearized independently, while its converter-voltage output is applied to the plant with an explicit one-time-step delay. This allows the controller matrices to be updated without changing the plant Norton matrix.
The PR includes:
DP::VTypeSplitSSNCompbase implementing packed-real DP states, trapezoidal discretization, controller relinearization, delayed controller output, and Norton-equivalent coupling;DP::Ph3::SSN_GFL_Split, a positive-sequence grid-following inverter with PLL, filtered-power control, current control, analytical controller Jacobians, and steady-state initialization;The example compares the extracted 32-state discrete model against an independently assembled delayed dq0 reference.