diff --git a/docs/user_manual/sphinx/Section10-Examples.rst b/docs/user_manual/sphinx/Section10-Examples.rst index 348d3e816..b51e41a97 100644 --- a/docs/user_manual/sphinx/Section10-Examples.rst +++ b/docs/user_manual/sphinx/Section10-Examples.rst @@ -1305,7 +1305,10 @@ An example contingency application has been included in the contingency analysis directory. This contingency analysis is simpler than the one available under the ``applications`` directory and provides a relatively compact demonstration of some of the advanced features of -GridPACK. This application is built entirely around the power flow +GridPACK. The production driver under ``src/applications/contingency_analysis`` +adds CSV outputs, monitor filters, and PSS/E-aligned options; its +runtime configuration is documented in the *Contingency Analysis Module* +section of Section 9. This application is built entirely around the power flow module, so it has no network component classes of its own. The main functionality is located in the ``CADriver`` class that consists of two methods (other than the constructor and destructor). One function is diff --git a/docs/user_manual/sphinx/Section9-ApplicationModules.rst b/docs/user_manual/sphinx/Section9-ApplicationModules.rst index d7fbe3645..3419d6024 100644 --- a/docs/user_manual/sphinx/Section9-ApplicationModules.rst +++ b/docs/user_manual/sphinx/Section9-ApplicationModules.rst @@ -135,6 +135,16 @@ adjusted by one discrete step per controller iteration, bounded by number of tap positions (NTP), or read directly from the STEP field when available. Cycle detection prevents tap hunting. +Remote voltage regulation (IREG) is enabled automatically for any +generator whose PSS/E ``IREG`` field points to a bus other than its +own. The controlled remote bus is treated as PV (held at the +generator's scheduled voltage ``VS``) and the local bus is solved +normally; the source bus voltage is adjusted across the controller loop +until the remote bus tracks setpoint. No XML option is required. When +the remote bus is the swing bus the regulation is dropped with a +warning, and when multiple generators on the same source bus include a +locally-regulating unit (``IREG=0``) local control wins. + The ``AreaInterchange`` parameter (default ``false``) enables area interchange control. When enabled, after the inner controller loop converges, the solver computes actual MW exports for each area by @@ -537,6 +547,255 @@ Again, this may be useful in contingency calculations where multiple calculations are run on the same network and it is desirable that they all start with the same initial condition. +Contingency Analysis Module +--------------------------- + +The production contingency analysis driver lives in +``src/applications/contingency_analysis``. It is built on +``PFAppModule`` and adds an MPI task manager, a contingency parser, +monitor filtering, and per-(contingency, branch) CSV output for +downstream statistical analysis. Configuration goes under a +``Contingency_analysis`` block in the input deck alongside the standard +``Powerflow`` block. The XML reader for the contingency list itself is +covered in Section 10; this section documents the runtime options. + +Input options +~~~~~~~~~~~~~ + ++----------------------------+--------+--------------------------------------------------------------+ +| Option | Default| Description | ++============================+========+==============================================================+ +| ``contingencyList`` | (none) | Path to the contingency XML. May be combined with | +| | | ``FullBranchN1`` / ``FullGeneratorN1`` for N-1 + custom N-K. | ++----------------------------+--------+--------------------------------------------------------------+ +| ``FullBranchN1`` | false | Auto-generate N-1 over every in-service branch. | ++----------------------------+--------+--------------------------------------------------------------+ +| ``FullGeneratorN1`` | false | Auto-generate N-1 over every in-service generator. | ++----------------------------+--------+--------------------------------------------------------------+ +| ``groupSize`` | 1 | Deprecated; ignored (forced to 1). An outaged branch can | +| | | straddle a multi-rank partition. Add ranks for parallelism. | ++----------------------------+--------+--------------------------------------------------------------+ +| ``minVoltage`` | 0.9 | Lower voltage limit (pu) for violation checks. | ++----------------------------+--------+--------------------------------------------------------------+ +| ``maxVoltage`` | 1.1 | Upper voltage limit (pu) for violation checks. | ++----------------------------+--------+--------------------------------------------------------------+ +| ``qlim`` | false | Enable PV→PQ Q-limit enforcement during the contingency | +| | | solve. Honors ``qlimDeadband`` from the ``Powerflow`` block. | ++----------------------------+--------+--------------------------------------------------------------+ +| ``printCalcFiles`` | true | Write per-contingency text output (``.out``). | ++----------------------------+--------+--------------------------------------------------------------+ +| ``outputFormat`` | text | ``text`` | ``json`` | ``csv`` | ``csv_flat`` | ``csv_delta``.| ++----------------------------+--------+--------------------------------------------------------------+ +| ``outputFile`` | | Base name (prefix) for the structured output files. | +| | results| | ++----------------------------+--------+--------------------------------------------------------------+ +| ``writeStats`` | true | Emit the StatBlock summary ``.txt`` files. Set ``false`` to | +| | | skip the per-case StatBlock work when CSV is sufficient. | ++----------------------------+--------+--------------------------------------------------------------+ +| ``contingencyRating`` | A | Loading% denominator across every CA output. ``A``, ``B``, | +| | | or ``C`` with A→B→C fallback. Default ``A`` matches PW / | +| | | PSS/E ACCC. ``base_rate_mva`` always uses rate-A. | ++----------------------------+--------+--------------------------------------------------------------+ +| ``monitorBranchesFile`` | (none) | Path to a CSV allowlist (``from_bus,to_bus,ckt`` rows). | +| | | When set, area/kV gates are ignored. | ++----------------------------+--------+--------------------------------------------------------------+ +| ``monitorAreas`` | (none) | Space-separated PSS/E area numbers. Branch passes if | +| | | **either endpoint** is in the set (catches tie-lines). | ++----------------------------+--------+--------------------------------------------------------------+ +| ``monitorKvMin`` | 0 | Lower kV threshold; branch passes if | +| | | ``max(kv_from, kv_to) >= monitorKvMin``. 0 disables. | ++----------------------------+--------+--------------------------------------------------------------+ +| ``monitorKvMax`` | 0 | Upper kV threshold; branch passes if | +| | | ``max(kv_from, kv_to) <= monitorKvMax``. 0 disables. | ++----------------------------+--------+--------------------------------------------------------------+ + +Violation / ranking options driving ``_violations.csv`` and +``_summary.json``: + +* ``violationSeverityThreshold`` (``1.0``) — loading% > threshold × 100 + flags a branch violation. +* ``topN`` (``20``, clamped ``[1,10000]``) — cap on ranked arrays. +* ``piBranchWeight`` / ``piVoltageWeight`` (both ``1.0``) — weights on + the branch-PI and voltage-PI terms in the composite PI ranking. + +Monitor filter precedence: When ``monitorBranchesFile`` is set, area/kV +options are ignored and a warning is logged. Otherwise ``monitorAreas`` +AND the kV bounds combine. With no filter set every branch is monitored. +Circuit IDs are matched with whitespace trimmed on both sides so PSS/E ckt +strings with leading or trailing spaces compare equal to the unpadded form. + +A complete annotated example is shipped at +``src/applications/data_sets/input/ca/input_14_filters_example.xml`` +with a sample monitor list ``monitor_branches_14.csv``. + +CSV outputs (``outputFormat=csv_flat`` / ``csv_delta``) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +When ``outputFormat`` is ``csv_flat`` or ``csv_delta`` the driver writes +per-(contingency, monitored branch) rows for downstream analysis instead +of the aggregated ``.txt`` files. All file names use the value of +``outputFile`` as a prefix; the default prefix is ``ca_results``. Both +formats also work for generator contingencies — column 3 (``type``) +labels what was tripped and column 31 (``cont_event_facility``) +identifies the tripped element. + +``_delta.csv`` (``csv_delta`` only) — wide form, one row +per (contingency, monitored branch) joining base and contingency state +on the same row. Most downstream consumers prefer this because each row +is self-contained. + ++-----+----------------------------+----------------------------------------------------------+ +| # | Column | Notes | ++=====+============================+==========================================================+ +| 1 | ``event_idx`` | 0 = base case, 1..N = contingencies in input-deck order. | ++-----+----------------------------+----------------------------------------------------------+ +| 2 | ``contingency`` | Contingency name (``base_case`` for the base row). | ++-----+----------------------------+----------------------------------------------------------+ +| 3 | ``type`` | ``branch`` or ``generator`` — what was tripped. | ++-----+----------------------------+----------------------------------------------------------+ +| 4–6 | ``from_bus``, ``to_bus``, | Identity of the **monitored branch** (not the tripped | +| | ``ckt`` | element). | ++-----+----------------------------+----------------------------------------------------------+ +| 7–8 | ``base_kv_from``, | Endpoint base kV. | +| | ``base_kv_to`` | | ++-----+----------------------------+----------------------------------------------------------+ +| 9–10| ``area_from``, | PSS/E area numbers. | +| | ``area_to`` | | ++-----+----------------------------+----------------------------------------------------------+ +| 11 | ``base_rate_mva`` | Always rate-A (PSS/E "normal" rating). | ++-----+----------------------------+----------------------------------------------------------+ +| 12 | ``cont_rate_mva`` | Rating selected by ``contingencyRating``. | ++-----+----------------------------+----------------------------------------------------------+ +| 13–14| ``base_p_mw``, | Real-power flow before / after. | +| | ``cont_p_mw`` | | ++-----+----------------------------+----------------------------------------------------------+ +| 15–16| ``base_q_mvar``, | Reactive-power flow before / after. | +| | ``cont_q_mvar`` | | ++-----+----------------------------+----------------------------------------------------------+ +| 17–18| ``base_mva``, | ``sqrt(P² + Q²)`` before / after. | +| | ``cont_mva`` | | ++-----+----------------------------+----------------------------------------------------------+ +| 19 | ``base_loading_pct`` | ``base_mva / base_rate_mva × 100``. | ++-----+----------------------------+----------------------------------------------------------+ +| 20 | ``cont_loading_pct`` | ``cont_mva / cont_rate_mva × 100``. | ++-----+----------------------------+----------------------------------------------------------+ +| 21–22| ``v_from_base``, | From-bus voltage magnitude (pu). | +| | ``v_from_cont`` | | ++-----+----------------------------+----------------------------------------------------------+ +| 23–24| ``v_to_base``, | To-bus voltage magnitude (pu). | +| | ``v_to_cont`` | | ++-----+----------------------------+----------------------------------------------------------+ +| 25–28| ``ang_from_base``, | Endpoint bus angles (deg). | +| | ``ang_from_cont``, | | +| | ``ang_to_base``, | | +| | ``ang_to_cont`` | | ++-----+----------------------------+----------------------------------------------------------+ +| 29–30| ``d_angle_base``, | ``ang_from − ang_to`` before / after. | +| | ``d_angle_cont`` | | ++-----+----------------------------+----------------------------------------------------------+ +| 31 | ``cont_event_facility`` | The tripped element (e.g. ``[area] from to ckt`` for a | +| | | branch trip, ``gen `` for a gen trip). | ++-----+----------------------------+----------------------------------------------------------+ + +``_flat.csv`` (``csv_flat`` only) — long form, one row per +(case, branch). Columns: ``event_idx, contingency, from_bus, to_bus, +ckt, p_from_mw, q_from_mvar, mva_from, rate_mva, loading_percent, viol, +v_from_pu, v_to_pu, ang_from_deg, ang_to_deg``. ``rate_mva`` is rate-A +on the base row and the configured ``contingencyRating`` on contingency +rows. + +``_buses.csv`` (both CSV formats) — bus metadata sidecar so +the per-branch files can stay narrow. Columns: ``bus_id, bus_name, +base_kv, area, zone, owner, area_name, zone_name, owner_name``. + +``_convergence.csv`` (every ``outputFormat``) — one row per +contingency. Columns: ``event_idx, contingency, type, converged, +iterations, final_tolerance, max_p_bus, max_p_mismatch, max_q_bus, +max_q_mismatch, status_code``. ``converged`` is ``true`` iff +``status_code == "OK"``; ``status_code`` is one of ``OK`` / ``ISLANDED`` +/ ``NO_SLACK`` / ``DIVERGED`` / ``SLACK_OVERLOAD``. Failed rows appear +here even when omitted from ``_delta.csv`` / ``_flat.csv``. + +``_violations.csv`` (every ``outputFormat``) — streamed row +per branch/voltage violation. Columns: ``event_idx, contingency, type, +element, mva_or_vpu, rate_or_limit, loading_percent, base_mva, delta, +severity``. Emitted when ``loading_percent > violationSeverityThreshold +× 100`` (branch) or ``v_pu`` outside ``[minVoltage, maxVoltage]`` +(voltage). ``loading_percent`` divides by the rate picked under +``contingencyRating``. ``severity = "critical"`` when branch loading +≥ 105 % or |Δv| ≥ 0.05 pu, else ``"warning"``. + +``_summary.json`` (every ``outputFormat``) — end-of-run +aggregate: ``total_contingencies``, ``converged``, ``diverged = +total − converged``, plus the per-status split ``islanded``, +``no_slack``, ``solver_diverged``, ``slack_overload`` (sum to +``diverged``). Also ``with_branch_violation`` / ``with_voltage_violation`` +/ ``violation_rows``; single-element extremes ``worst_loading`` / +``worst_voltage_low`` / ``worst_voltage_high`` (``null`` if none); +name arrays ``contingencies_with_branch_violation`` and +``contingencies_with_voltage_violation``; a length-``topN`` ranked +``top_severe_contingencies`` list (violated cases first by severity, +non-violated after by composite PI); and an echo of the run config +(``contingency_rating``, ``voltage_limit_low``/``high``, +``severity_threshold``). + +When monitor filters are active the data-row count of ``_delta.csv`` / +``_flat.csv`` equals ``|monitored branches| × |converged +contingencies|``. + +Output ordering +~~~~~~~~~~~~~~~ + +Rows are not sorted by contingency. The driver streams each MPI rank's +results to a per-rank ``.part`` file and rank 0 concatenates them in +rank order, so the final file is grouped by rank and ordered by +completion within each rank. Column 1 (``event_idx``) preserves +input-deck order — sort downstream if needed:: + + ( head -1 my_run_delta.csv && \ + tail -n +2 my_run_delta.csv | sort -t, -k1,1n ) > my_run_delta.sorted.csv + +Aggregated ``.txt`` outputs (``writeStats=true``) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +When ``writeStats=true`` the driver also writes the StatBlock summary +files: ``vmag.txt``, ``vmag_mm.txt``, ``vang.txt``, ``vang_mm.txt``, +``pgen.txt``, ``pgen_mm.txt``, ``qgen.txt``, ``qgen_mm.txt``, +``pflow.txt``, ``pflow_mm.txt``, ``qflow.txt``, ``qflow_mm.txt``, +``perf_mm.txt``, ``perf_sum.txt``, ``line_flt_cnt.txt``. The file +``pq_change_cnt.txt`` is added when ``qlim=true``. These hold per-element +statistics (mean / RMS / min / max) **across all contingencies**, not +per-contingency rows; their layout is described in the application +``README.md``. Set ``writeStats=false`` to skip them when CSV output is +sufficient. + +Contingency XML aliases +~~~~~~~~~~~~~~~~~~~~~~~ + +The element names in a contingency XML accept PSS/E-aligned aliases: + ++----------------------------+------------+-------------------------------+ +| Element | Alias | Holds | ++============================+============+===============================+ +| ``contingencyLineNames`` | ``CKT`` | Branch circuit ID (PSS/E | +| | | ``CKT`` field). | ++----------------------------+------------+-------------------------------+ +| ``contingencyGenerators`` | ``GenID`` | Generator ID (PSS/E ``ID`` | +| | | field). | ++----------------------------+------------+-------------------------------+ + +Either name parses; mix-and-match within a file is fine. + +Advanced behavior +~~~~~~~~~~~~~~~~~ + +The driver also includes automatic slack-bus transfer when the slack +generator is the tripped element, a slack-capacity check after each +solve, and island / lone-bus detection when a branch trip splits the +network. The underlying mechanisms (``checkAndTransferSlack``, +``restoreSlack``, ``checkSlackCapacity``, ``getIslandCount``, +``hasLoneBus``) are exposed by ``PFAppModule`` and described above. + State Estimation Module ----------------------- diff --git a/src/applications/components/pf_matrix/pf_components.cpp b/src/applications/components/pf_matrix/pf_components.cpp index df8efd8be..05e0026fa 100755 --- a/src/applications/components/pf_matrix/pf_components.cpp +++ b/src/applications/components/pf_matrix/pf_components.cpp @@ -69,6 +69,7 @@ // Static member initialization gridpack::powerflow::InitStartMode gridpack::powerflow::PFBus::p_initStartMode = INIT_START_WARM; bool gridpack::powerflow::PFBus::p_qlim = true; +double gridpack::powerflow::PFBus::p_qlim_deadband = 0.1; std::vector gridpack::powerflow::PFBus::p_qlimWarnings; /** @@ -87,6 +88,11 @@ void gridpack::powerflow::PFBus::setQlim(bool qlim) p_qlim = qlim; } +void gridpack::powerflow::PFBus::setQlimDeadband(double db) +{ + p_qlim_deadband = db; +} + /** * Clear accumulated Q limit warning messages */ @@ -453,7 +459,7 @@ bool gridpack::powerflow::PFBus::chkQlim(double q_deadband) // Check if Q requirement can be met. // q_deadband avoids switching buses that are only marginally over their Q // limit due to floating-point differences. Configurable via XML qlimDeadband - // (default 0.1 Mvar, matching PW's 0.1 MVA convergence tolerance). + // (default 0.1 Mvar) bool need_pv_to_pq = false; char warnBuf[256]; if (Q_required > Q_max_total + q_deadband) { @@ -1077,11 +1083,7 @@ void gridpack::powerflow::PFBus::load( } } } - // Warm-start Q-limit pre-saturation: if scheduled QG is already at QMAX or QMIN - // (within 0.1 Mvar), start this bus as PQ immediately. This matches PW behavior - // where generators already at their limits in the input data are treated as PQ, - // preventing IREG from driving them to impossible Q requirements. - // Only applies for warm start with qlim enabled. + // Warm-start Q-limit pre-saturation: scheduled QG at QMAX/QMIN -> start as PQ. if (p_isPV && p_qlim && p_initStartMode != INIT_START_FLAT) { double total_qg = 0.0, total_qmax = 0.0, total_qmin = 0.0; for (i = 0; i < p_ngen; i++) { @@ -1091,9 +1093,23 @@ void gridpack::powerflow::PFBus::load( total_qmin += p_qmin[i]; } } - const double Q_init_tol = 0.1; // Mvar — matches PW convergence tolerance - if (total_qg >= total_qmax - Q_init_tol || total_qg <= total_qmin + Q_init_tol) { + if (total_qg >= total_qmax - p_qlim_deadband + || total_qg <= total_qmin + p_qlim_deadband) { p_isPV = false; + // Pick the saturated limit from V vs VS (V>VS => QMIN, V QMAX), not + // scheduled QG which can be inconsistent. Local regulation only. + if (p_ireg_remote_bus == 0) { + double vset = 0.0; int nv = 0; + for (i = 0; i < p_ngen; i++) + if (p_gstatus[i] == 1) { vset += p_vs[i]; nv++; } + if (nv > 0) vset /= nv; + const double V_init_tol = 1.0e-4; + if (p_voltage > vset + V_init_tol) { + for (i = 0; i < p_ngen; i++) if (p_gstatus[i] == 1) p_qg[i] = p_qmin[i]; + } else if (p_voltage < vset - V_init_tol) { + for (i = 0; i < p_ngen; i++) if (p_gstatus[i] == 1) p_qg[i] = p_qmax[i]; + } + } } } @@ -1821,6 +1837,8 @@ bool gridpack::powerflow::PFBus::serialWrite(char *string, const int bufsize, sprintf(string, "%6d %20.12e %20.12e %d %d\n", getOriginalIndex(),0.0,0.0,use_vmag,changed); } else if (!strcmp(signal,"ca")) { + // Match checkVoltageViolations(): skip ignored buses. + if (p_ignore) return false; double pi = 4.0*atan(1.0); double angle = p_a*180.0/pi; bool found = false; @@ -2495,6 +2513,39 @@ int gridpack::powerflow::PFBus::getZone() return p_zone; } +/** + * Get owner number for bus + * @return bus owner number (0 if not set) + */ +int gridpack::powerflow::PFBus::getOwner() +{ + int owner = 0; + if (p_data) p_data->getValue(BUS_OWNER, &owner); + return owner; +} + +/** + * Get base voltage for bus in kV + * @return base kV (0.0 if not set) + */ +double gridpack::powerflow::PFBus::getBaseKV() +{ + double basekv = 0.0; + if (p_data) p_data->getValue(BUS_BASEKV, &basekv); + return basekv; +} + +/** + * Get bus name string + * @return bus name (empty if not set) + */ +std::string gridpack::powerflow::PFBus::getBusName() +{ + std::string name; + if (p_data) p_data->getValue(BUS_NAME, &name); + return name; +} + /** * Evaluate diagonal block of Jacobian for power flow calculation and * return result as an array of real values @@ -3522,7 +3573,16 @@ gridpack::ComplexType gridpack::powerflow::PFBranch::getComplexPower( { gridpack::ComplexType vi, vj, Yii, Yij, s; s = ComplexType(0.0,0.0); - gridpack::powerflow::PFBus *bus1 = + // Out-of-service line: getLineElements does not gate on status, so + // return zero to avoid a phantom flow from stale admittance. + int bsize = p_branch_status.size(); + for (int i=0; i(getBus1().get()); vi = bus1->getComplexVoltage(); gridpack::powerflow::PFBus *bus2 = @@ -3543,6 +3603,14 @@ gridpack::ComplexType gridpack::powerflow::PFBranch::getReversePower( { gridpack::ComplexType vi, vj, Yjj, Yji, s; s = ComplexType(0.0,0.0); + // Out-of-service line: see getComplexPower. + int bsize = p_branch_status.size(); + for (int i=0; i(getBus1().get()); vi = bus1->getComplexVoltage(); @@ -3562,10 +3630,42 @@ gridpack::ComplexType gridpack::powerflow::PFBranch::getReversePower( * routine what about kind of information to write * @return true if branch is contributing string to output, false otherwise */ +std::string gridpack::powerflow::PFBranch::s_contingencyRating = "A"; + +void gridpack::powerflow::PFBranch::setContingencyRating( + const std::string& rating) +{ + if (rating == "A" || rating == "B" || rating == "C") { + s_contingencyRating = rating; + } else { + s_contingencyRating = "A"; + } +} + +std::string gridpack::powerflow::PFBranch::getContingencyRating() +{ + return s_contingencyRating; +} + +// A->B->C fallback when the picked tier is zero. +double gridpack::powerflow::PFBranch::pickBranchRating(int elemIdx) const +{ + if (elemIdx < 0 || elemIdx >= static_cast(p_rateA.size())) return 0.0; + double a = p_rateA[elemIdx]; + double b = (elemIdx < static_cast(p_rateB.size())) ? p_rateB[elemIdx] : 0.0; + double c = (elemIdx < static_cast(p_rateC.size())) ? p_rateC[elemIdx] : 0.0; + if (s_contingencyRating == "A") return a; + if (s_contingencyRating == "B") return (b > 0.0) ? b : a; + // "C" + if (c > 0.0) return c; + if (b > 0.0) return b; + return a; +} + bool gridpack::powerflow::PFBranch::serialWrite(char *string, const int bufsize, const char *signal) { - char buf[128]; + char buf[256]; gridpack::powerflow::PFBus *bus1 = dynamic_cast(getBus1().get()); gridpack::powerflow::PFBus *bus2 @@ -3636,18 +3736,33 @@ bool gridpack::powerflow::PFBranch::serialWrite(char *string, const int bufsize, bool found = false; int ilen = 0; for (i=0; i(p_ignore.size()) && p_ignore[i]) continue; s = getComplexPower(tags[i]); double p = real(s); double q = imag(s); + gridpack::ComplexType s2 = getReversePower(tags[i]); + double p2 = real(s2); + double q2 = imag(s2); if (!p_branch_status[i]) p = 0.0; if (!p_branch_status[i]) q = 0.0; + if (!p_branch_status[i]) p2 = 0.0; + if (!p_branch_status[i]) q2 = 0.0; if (bus1->isIsolated() || bus2->isIsolated()) p=0.0; if (bus1->isIsolated() || bus2->isIsolated()) q=0.0; - double S = sqrt(p*p+q*q); - if (S > p_rateA[i] && p_rateA[i] != 0.0){ - sprintf(buf, " %6d %6d %s %12.6f %12.6f %8.2f %8.2f%s\n", + if (bus1->isIsolated() || bus2->isIsolated()) p2=0.0; + if (bus1->isIsolated() || bus2->isIsolated()) q2=0.0; + double Sfrom = sqrt(p*p+q*q); + double Sto = sqrt(p2*p2+q2*q2); + // One loading value, measured at the more heavily loaded end, as + // PowerWorld reports it and as _violations.csv/_branches.csv compute it. + double S = (Sfrom > Sto) ? Sfrom : Sto; + double rate = pickBranchRating(i); + if (S > rate && rate != 0.0){ + sprintf(buf, "%10d%10d%6s%13.6f%13.6f%13.6f%13.6f%13.6f%13.6f" + "%11.2f%10.2f%s\n", getBus1OriginalIndex(),getBus2OriginalIndex(),tags[i].c_str(), - p,q,p_rateA[i],S/p_rateA[i]*100,"%"); + p,q,Sfrom,p2,q2,Sto,rate,S/rate*100,"%"); int len = strlen(buf); if (ilen + len < bufsize) { sprintf(string,"%s",buf); diff --git a/src/applications/components/pf_matrix/pf_components.hpp b/src/applications/components/pf_matrix/pf_components.hpp index 72d250b8c..08a41cf39 100644 --- a/src/applications/components/pf_matrix/pf_components.hpp +++ b/src/applications/components/pf_matrix/pf_components.hpp @@ -481,6 +481,24 @@ class PFBus */ int getZone(); + /** + * Get owner number for bus + * @return bus owner number (0 if not set) + */ + int getOwner(); + + /** + * Get base voltage for bus in kV + * @return base kV (0.0 if not set) + */ + double getBaseKV(); + + /** + * Get bus name string + * @return bus name (empty if not set) + */ + std::string getBusName(); + /** * Evaluate diagonal block of Jacobian for power flow calculation and return * result as an array of real values @@ -610,6 +628,7 @@ class PFBus */ static void setInitStartMode(InitStartMode mode); static void setQlim(bool qlim); + static void setQlimDeadband(double db); /** * Clear accumulated Q limit warning messages @@ -691,6 +710,7 @@ class PFBus static std::vector p_qlimWarnings; static InitStartMode p_initStartMode; static bool p_qlim; + static double p_qlim_deadband; double p_shunt_gs; double p_shunt_bs; bool p_shunt; @@ -970,6 +990,12 @@ class PFBranch */ double getBranchRatingC(std::string tag); + // Contingency rating tier ("A"|"B"|"C") used by serialWrite("flow",...). + // Default "A" preserves non-CA behavior; A->B->C fallback if picked=0. + static void setContingencyRating(const std::string& rating); + static std::string getContingencyRating(); + double pickBranchRating(int elemIdx) const; + /** * Get list of line IDs * @return list of line identifiers @@ -1057,6 +1083,8 @@ class PFBranch // LTC (Load Tap Changer) control variables bool p_hasLTC; // true if this branch has an LTC-controlled transformer int p_ltc_elem; // index of the LTC element within this branch + static std::string s_contingencyRating; + int p_ltc_code; // control mode (1=voltage, 0=off) int p_ltc_cont; // controlled bus number bool p_ltc_cont_is_to; // true if controlled bus is the to-bus (tap direction reversal) diff --git a/src/applications/contingency_analysis/README.md b/src/applications/contingency_analysis/README.md index f95934f29..61adbd360 100644 --- a/src/applications/contingency_analysis/README.md +++ b/src/applications/contingency_analysis/README.md @@ -35,16 +35,110 @@ When combined, duplicates from the file are automatically skipped. | Option | Description | Default | |--------|-------------|---------| -| `groupSize` | Number of MPI processes per contingency (parallelization) | 1 | +| `groupSize` | Deprecated; ignored. Each contingency runs on one rank (a group larger than one lets an outaged branch straddle the partition and fails the solve). Add MPI ranks to run more contingencies at once | 1 (forced) | | `printCalcFiles` | Write detailed output for each contingency | true | | `minVoltage` | Minimum voltage threshold for violations (p.u.) | 0.9 | | `maxVoltage` | Maximum voltage threshold for violations (p.u.) | 1.1 | | `qlim` | Enable reactive power limit enforcement (PV to PQ bus conversion) | false | +| `outputFormat` | `text` / `json` / `csv` / `csv_flat` / `csv_delta` | `text` | +| `outputFile` | Base name for output files | `ca_results` | +| `writeStats` | Emit StatBlock summary files (vmag.txt etc.). Set false to skip and avoid the per-case StatBlock work | true | +| `contingencyRating` | Loading% denominator across every CA output: `A`, `B`, or `C` with A→B→C fallback. Default `A` matches PW / PSS/E ACCC. `base_rate_mva` always uses rate-A | `A` | +| `monitorBranchesFile` | Path to a CSV allowlist (`from_bus,to_bus,ckt`). When set, overrides the area/kV gates (TARA / PSS/E convention). Applies to every output format | (unset) | +| `monitorAreas` | Space-separated list of PSS/E area numbers. A branch is monitored if **either endpoint** is in the set; a bus if its own area is. Applies to every output format | (unset) | +| `monitorKvMin` | Lower kV threshold; branch passes if `max(kv_from, kv_to) >= monitorKvMin`, bus if its base kV does. Applies to every output format | 0 (unbounded) | +| `monitorKvMax` | Upper kV threshold; branch passes if `max(kv_from, kv_to) <= monitorKvMax`, bus if its base kV does. Applies to every output format | 0 (unbounded) | + +### Monitor filters (all output formats) + +The monitor filters select which **buses and branch elements are reported**; +they never change which contingencies are simulated or how the network is +solved. All filters are optional — unset means "monitor everything". +`monitorBranchesFile` is authoritative when set; otherwise `monitorAreas` and +the kV bounds AND together. + +The same monitored set is used by every output, regardless of `outputFormat`: + +| Output | What the filter does | +|---|---| +| `.out` per-contingency text files (`printCalcFiles=true`) | Bus / branch violation listings and the "No violation" verdict cover monitored elements only | +| `_violations.csv` | Branch and voltage rows for monitored elements only | +| `json` / `csv` (`_buses.csv`, `_branches.csv`, `_generators.csv`, `.json`) | Base-case and contingency rows for monitored buses, branches and generators on monitored buses | +| `csv_flat` / `csv_delta` | One row per (contingency, monitored branch) | +| StatBlock `.txt` files (`vmag.txt`, `pflow.txt`, ...) | Rows for monitored buses / branches / generators only | +| `_summary.json` | Violation counters, worst-of values, performance indices and rosters accrue over monitored elements only | +| `_buses.csv` metadata sidecar (csv_flat / csv_delta), `_contingencies.csv`, `_convergence.csv` | Not filtered: these are lookup / bookkeeping tables and stay complete | + +The bus rule mirrors the branch rule: with an allowlist a bus is monitored when +it is an endpoint of an allowlisted branch; with area / kV gates it is monitored +when its own area is in `monitorAreas` and its own base kV lies inside the kV +bounds. Non-monitored elements are flagged "ignore" on the network after the +base-case solve, which is the same mechanism the driver already uses to exclude +buses that violate limits in the base case from the contingency checks. + +At startup the driver prints the effective filter and a count such as +`Monitor filter: 20 of 179 branch elements and 10 of 118 buses monitored`. +A filter that matches nothing produces a warning and empty outputs; check the +area numbers and kV levels against the case (legacy v23 RAW files often carry +area 1 and base kV 0 for every bus, so area / kV gates cannot select anything +there). + +```xml + + + +monitor_branches.csv + + +11 12 19 +100.0 +500.0 +``` + +`monitorBranchesFile` is a CSV of `from_bus,to_bus,ckt` rows (header +optional; `#` is a line comment). When set, area/kV options are ignored +and a warning is logged. + +`monitorAreas` matches branches with either endpoint in the set (catches +tie-lines). `monitorKvMin/Max` gate on `max(kv_from, kv_to)` so a 138/13.8 +step-down counts as 138. + +`contingencyRating` (`A` | `B` | `C`, default `A`) sets the loading% +denominator for every CA output. Default `A` matches PowerWorld / PSS/E +ACCC. `base_rate_mva` always uses rate-A. A→B→C fallback if the picked +tier is zero/missing. + +A complete annotated example is in +`src/applications/data_sets/input/ca/input_14_filters_example.xml` with a +sample monitor file `monitor_branches_14.csv` in the same directory. + +### Output ordering + +Rows are not sorted by contingency. The driver distributes contingencies +across MPI ranks and streams each rank's results to its own `.part` file; +rank 0 concatenates in rank order, so the final file is grouped by rank +and ordered by completion within each rank. Column 1 (`event_idx`) +preserves input-deck order — sort downstream if needed: + +```bash +( head -1 my_run_delta.csv && tail -n +2 my_run_delta.csv | sort -t, -k1,1n ) > my_run_delta.sorted.csv +``` + ### Contingency File Format See `contingencies_nk_example.xml` for examples of N-1, N-2, and N-3 contingency definitions. +The line-tag and generator-id elements accept PSS/E-aligned aliases for clarity: + +| Element | Alias | Holds | +|---|---|---| +| `` | `` | Branch circuit ID (PSS/E `CKT` field) | +| `` | `` | Generator ID (PSS/E `ID` field) | + +Either name works; mix-and-match within the same file is fine. Both legacy +files and new files using the PSS/E names continue to parse without changes. + --- ## Advanced Features @@ -64,13 +158,6 @@ After the contingency analysis completes, the slack bus is restored to its origi After the power flow solves, the application checks if the slack bus generator output exceeds its Pmax rating. If the required generation exceeds capacity, the contingency is marked as failed with a warning message: -``` -WARNING: Slack bus 80 generator output (475.3 MW) exceeds capacity (400.0 MW) -Insufficient generation capacity for contingency GN_69_1 -``` - -This ensures realistic results - a contingency that requires more generation than available capacity is properly flagged as a failure. - ### Island Detection The application detects network islands (disconnected portions) caused by branch contingencies: @@ -90,23 +177,93 @@ for contingencies that ran to completion are included. Calculations that failed either because of a numerical instability or because the calculations failed to converge are not included in the results. The output files are described below. -**success.txt**: This file summarizes the results of each contingency and -reports 1) whether the contingency calculation successfully ran to completion, -2) whether a violation was found (bus, branch, or both), and 3) whether any -buses were isolated. +### CSV outputs (`outputFormat=csv_flat` / `csv_delta`) + +When `outputFormat` is set to `csv_flat` or `csv_delta`, the application writes +per-(contingency, branch) rows for downstream statistical analysis instead of +the aggregated `.txt` files described later in this section. All file names +below use the value of `outputFile` as a prefix; the default prefix is +`ca_results`. + +**`_delta.csv`** *(`outputFormat=csv_delta` only)* — wide-form, +one row per (contingency, monitored branch) joining base + contingency state +on the same row. This is the format most downstream consumers prefer because +each row is self-contained (no separate base-case join needed). + +| # | Column | Notes | +|---|---|---| +| 1 | `event_idx` | 0 = base case (only if a base row is emitted), 1..N = contingencies in input-deck order | +| 2 | `contingency` | contingency name from the input XML (`base_case` for the base) | +| 3 | `type` | `branch` or `generator` — what kind of contingency was tripped | +| 4–6 | `from_bus`, `to_bus`, `ckt` | Identity of the **monitored branch** in the row (not the tripped element) | +| 7–8 | `base_kv_from`, `base_kv_to` | Endpoint base kV | +| 9–10 | `area_from`, `area_to` | PSS/E area numbers | +| 11 | `base_rate_mva` | Always rate-A (PSS/E "normal" rating) | +| 12 | `cont_rate_mva` | Rating selected by `contingencyRating` (default A, with A→B→C fallback if zero/missing) | +| 13–14 | `base_p_mw`, `cont_p_mw` | Real-power flow before / after contingency | +| 15–16 | `base_q_mvar`, `cont_q_mvar` | Reactive-power flow before / after | +| 17–18 | `base_mva`, `cont_mva` | `sqrt(P² + Q²)` before / after | +| 19 | `base_loading_pct` | `base_mva / base_rate_mva × 100` | +| 20 | `cont_loading_pct` | `cont_mva / cont_rate_mva × 100` | +| 21–22 | `v_from_base`, `v_from_cont` | From-bus voltage magnitude (pu) before / after | +| 23–24 | `v_to_base`, `v_to_cont` | To-bus voltage magnitude (pu) before / after | +| 25–28 | `ang_from_base`, `ang_from_cont`, `ang_to_base`, `ang_to_cont` | Bus angles (deg) | +| 29–30 | `d_angle_base`, `d_angle_cont` | `ang_from − ang_to` before / after | +| 31 | `cont_event_facility` | Identifier of the **tripped element** in this contingency (e.g. `[area] from to ckt` for branch trips, `gen ` for gen trips) | + +**`_flat.csv`** *(`outputFormat=csv_flat` only)* — long-form, +one row per (case, branch). Columns: `event_idx, contingency, from_bus, +to_bus, ckt, p_from_mw, q_from_mvar, mva_from, rate_mva, loading_percent, +viol, v_from_pu, v_to_pu, ang_from_deg, ang_to_deg`. `rate_mva` is rate-A +on `event_idx=0` rows, the configured `contingencyRating` on contingency +rows. + +**`_buses.csv`** *(both `csv_flat` and `csv_delta`)* — bus +metadata sidecar so the per-branch files can stay narrow. Columns: +`bus_id, bus_name, base_kv, area, zone, owner, area_name, zone_name, +owner_name`. + +**`_convergence.csv`** *(every `outputFormat`)* — one row per +contingency. Columns: `event_idx, contingency, type, converged, +iterations, final_tolerance, max_p_bus, max_p_mismatch, max_q_bus, +max_q_mismatch, status_code`. `converged` is `true` iff `status_code == +"OK"`; `status_code` is one of `OK` / `ISLANDED` / `NO_SLACK` / +`DIVERGED` / `SLACK_OVERLOAD`. Failed rows appear here even though +they're omitted from `_delta.csv` / `_flat.csv`. + +When monitor filters are active, the data-row count of `_delta.csv` / +`_flat.csv` equals `|monitored branches| × |converged contingencies|`. + +### `_contingencies.csv` sidecar (every `outputFormat`) + +Lookup table decoding the `event_idx` used by every other output. Columns: +`event_idx, contingency, type, n_elements, from_bus, to_bus, circuit_id, +gen_bus, gen_id`. One row per contingency; `event_idx=0` is the base case +(`n_elements=0`, id columns blank). An N-1 row carries the single outaged +element in the id columns, so it joins directly against the other CSVs. A +multi-element (N-k) event stays on one row with `n_elements=k` and the +element ids `;`-separated inside the same columns, e.g. `from_bus=1;3`, +`to_bus=2;4`, `circuit_id=1;1`. Generator events fill `gen_bus`/`gen_id` +and leave the branch columns blank. -Example output: -``` -contingency: 1 success: true violation: none -contingency: 2 success: true violation: branch -contingency: 3 success: true violation: none warning: isolated -contingency: 4 success: false -``` +--- -- `success: true` - Power flow converged and slack capacity is within limits -- `success: false` - Power flow failed, island detected, or slack capacity exceeded -- `violation: none/bus/branch` - Whether voltage or thermal limits were violated -- `warning: isolated` - One or more buses were isolated (lone bus or island) +### Aggregated `.txt` outputs (`writeStats=true`, default) + +The remaining files in this section are produced by the StatBlock summary +pipeline, controlled by the `writeStats` option (default `true`). They +contain per-element statistics (mean / RMS / min / max) aggregated **across +all contingencies**, not per-contingency rows. Set `writeStats=false` to +skip them when csv_flat / csv_delta output is sufficient. + +The set of files emitted is fixed; their names are not configurable. With +`writeStats=true` you get: `vmag.txt`, `vmag_mm.txt`, `vang.txt`, +`vang_mm.txt`, `pgen.txt`, `pgen_mm.txt`, `qgen.txt`, `qgen_mm.txt`, +`pflow.txt`, `pflow_mm.txt`, `qflow.txt`, `qflow_mm.txt`, `perf_mm.txt`, +`perf_sum.txt`, `line_flt_cnt.txt`. The file `pq_change_cnt.txt` is also +written when `qlim=true`. Per-contingency convergence/status is reported +via the `_convergence.csv` sidecar described above, which is written for +every `outputFormat`. **vmag.txt**: This file contains the average value of the voltage magnitude for non-PV buses. It also contains the RMS fluctuations of the voltage magnitude diff --git a/src/applications/contingency_analysis/ca_driver.cpp b/src/applications/contingency_analysis/ca_driver.cpp index 88d385201..c38349db8 100644 --- a/src/applications/contingency_analysis/ca_driver.cpp +++ b/src/applications/contingency_analysis/ca_driver.cpp @@ -15,6 +15,11 @@ * - Q-limit support integration * @date 2026-01-31 * + * @updated Yousu Chen + * - csv_flat / csv_delta per-(contingency,branch) outputs + * - monitorBranchesFile / monitorAreas / monitorKvMin/Max filters (all formats) + * @date 2026-06-21 + * * @brief Driver for contingency analysis calculation that make use of the * powerflow module to implement individual power flow simulations for * each contingency. The different contingencies are distributed across @@ -29,10 +34,21 @@ #include "gridpack/utilities/results_exporter.hpp" #include "ca_driver.hpp" +#include #include - -#define USE_SUCCESS -#define USE_STATBLOCK +#include +#include +#include +#include +#include +#include +#include +#include + +// Statistical-summary output (vmag.txt, pflow.txt, etc.) used to be controlled +// by a USE_STATBLOCK build-time macro; it is now a runtime XML option, +// `Configuration.Contingency_analysis.writeStats`, defaulting to true to +// preserve existing behavior. // Sets up multiple communicators so that individual contingency calculations // can be run concurrently @@ -79,7 +95,9 @@ std::vector std::string buses; contingencies[idx]->get("contingencyLineBuses",&buses); std::string names; - contingencies[idx]->get("contingencyLineNames",&names); + if (!contingencies[idx]->get("CKT",&names)) { + contingencies[idx]->get("contingencyLineNames",&names); + } // Tokenize bus string to get a list of individual buses std::vector string_vec = utils.blankTokenizer(buses); // Convert buses from character strings to ints @@ -115,7 +133,9 @@ std::vector std::string buses; contingencies[idx]->get("contingencyBuses",&buses); std::string gens; - contingencies[idx]->get("contingencyGenerators",&gens); + if (!contingencies[idx]->get("GenID",&gens)) { + contingencies[idx]->get("contingencyGenerators",&gens); + } // Tokenize bus string to get a list of individual buses std::vector string_vec = utils.blankTokenizer(buses); std::vector bus_ids; @@ -199,7 +219,9 @@ std::vector } } - printf("Auto-generated %d N-1 branch contingencies\n", branch_count); + if (gridpack::parallel::Communicator().rank() == 0) { + printf("Auto-generated %d N-1 branch contingencies\n", branch_count); + } } // Generate N-1 generator contingencies @@ -234,7 +256,9 @@ std::vector } } - printf("Auto-generated %d N-1 generator contingencies\n", gen_count); + if (gridpack::parallel::Communicator().rank() == 0) { + printf("Auto-generated %d N-1 generator contingencies\n", gen_count); + } } return ret; @@ -366,9 +390,27 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) print_calcs = true; } } + // Statistical-summary output (vmag.txt, pflow.txt, etc. via StatBlock). + // Default true to preserve existing behavior; set false to skip the + // per-case StatBlock work and the 13 post-loop global writes. + bool write_stats = true; + if (cursor->get("writeStats",&tmp_bool)) { + util.toLower(tmp_bool); + write_stats = (tmp_bool != "false"); + } + // groupSize is forced to 1: an outaged branch may straddle a multi-rank + // partition. Scale by adding ranks, not by widening a group. if (!cursor->get("groupSize",&grp_size)) { grp_size = 1; } + if (grp_size != 1) { + if (world.rank() == 0) { + printf("WARNING: groupSize=%d is not supported (a contingency branch " + "may span the partition boundary); using groupSize=1\n", + grp_size); + } + grp_size = 1; + } if (!cursor->get("minVoltage",&Vmin)) { Vmin = 0.9; } @@ -377,29 +419,84 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) } // Check for Q limit violations (qlim: true=enabled, false=disabled) bool check_Qlim = cursor->get("qlim", true); - // Output format: "json", "csv", or "text" (default) + double qlim_deadband = cursor->get("qlimDeadband", 0.1); + // Output format: "text" (default), "json", "csv", "csv_flat", "csv_delta". std::string outputFormat = "text"; cursor->get("outputFormat", &outputFormat); + if (outputFormat != "text" && outputFormat != "json" && + outputFormat != "csv" && outputFormat != "csv_flat" && + outputFormat != "csv_delta") { + if (world.rank() == 0) { + printf("ERROR: unrecognized outputFormat='%s'. " + "Must be one of: text, json, csv, csv_flat, csv_delta. Aborting.\n", + outputFormat.c_str()); + } + world.barrier(); + MPI_Abort(static_cast(world), 1); + } std::string outputFile = "ca_results"; cursor->get("outputFile", &outputFile); + // Optional CSV allowlist (from_bus,to_bus,ckt). Empty -> emit all. + std::string monitorBranchesFile; + cursor->get("monitorBranchesFile", &monitorBranchesFile); + // Optional area/kV gates. Empty/zero/missing -> no restriction on that + // dimension. Filters AND together with monitorBranchesFile. + std::string monitorAreasStr; + cursor->get("monitorAreas", &monitorAreasStr); + double monitorKvMin = 0.0; + cursor->get("monitorKvMin", &monitorKvMin); + double monitorKvMax = 0.0; + cursor->get("monitorKvMax", &monitorKvMax); + std::set monitorAreas; + { + std::vector tok = util.blankTokenizer(monitorAreasStr); + for (size_t i = 0; i < tok.size(); i++) { + if (!tok[i].empty()) monitorAreas.insert(atoi(tok[i].c_str())); + } + } + // Any monitor filter configured; gates every output format. + bool haveMonitorFilter = !monitorAreas.empty() || + monitorKvMin > 0.0 || monitorKvMax > 0.0 || + !monitorBranchesFile.empty(); + // Loading% denominator for all CA outputs (.out, _violations.csv, JSON, + // csv_flat, csv_delta). A|B|C, default A to match PW/PSSE convention. + // A->B->C fallback if the requested tier is zero/missing. + std::string contingencyRating = "A"; + cursor->get("contingencyRating", &contingencyRating); + util.toUpper(contingencyRating); + if (contingencyRating != "A" && contingencyRating != "B" && + contingencyRating != "C") { + if (world.rank() == 0) { + printf("WARNING: contingencyRating='%s' not A/B/C; defaulting to A\n", + contingencyRating.c_str()); + } + contingencyRating = "A"; + } + // Severity threshold for violation reporting; loading% > threshold*100 + // is flagged. Default 1.0 (100% of rate). + double violationSeverityThreshold = 1.0; + cursor->get("violationSeverityThreshold", &violationSeverityThreshold); + if (violationSeverityThreshold <= 0.0) violationSeverityThreshold = 1.0; + // Cap on top_severe_contingencies and roster arrays in _summary.json. + int topN = 10; + cursor->get("topN", &topN); + if (topN < 1) topN = 1; + if (topN > 10000) topN = 10000; + // Weights for composite_pi = piBranchWeight*branch_pi + piVoltageWeight*voltage_pi. + double piBranchWeight = 1.0, piVoltageWeight = 1.0; + cursor->get("piBranchWeight", &piBranchWeight); + cursor->get("piVoltageWeight", &piVoltageWeight); + if (piBranchWeight < 0.0) piBranchWeight = 0.0; + if (piVoltageWeight < 0.0) piVoltageWeight = 0.0; + // Set static flag for PFBus class BEFORE network creation. // This controls how Q values are reported in output functions: // - When check_Qlim = false: output uses calculated Q from p_Qinj // - When check_Qlim = true: output uses p_qg (set by chkQlim()) gridpack::powerflow::PFBus::setQlim(check_Qlim); + gridpack::powerflow::PFBus::setQlimDeadband(qlim_deadband); gridpack::parallel::Communicator task_comm = world.divide(grp_size); - // Keep track of failed calculations -#ifdef USE_SUCCESS - std::vector contingency_idx; - std::vector contingency_success; - gridpack::parallel::GlobalVector ca_success(world); - std::vector contingency_violation; - gridpack::parallel::GlobalVector ca_violation(world); - std::vector contingency_isolated; - gridpack::parallel::GlobalVector ca_isolated(world); -#endif - // Create powerflow applications on each task communicator boost::shared_ptr pf_network(new gridpack::powerflow::PFNetwork(task_comm)); @@ -410,23 +507,976 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) pf_app.readNetwork(pf_network,config); // Finish initializing the network pf_app.initialize(); + + // Build (number -> name) lookup tables for area, zone, owner. + // Keyed on the PSS/E-assigned number (not contiguous), used when + // emitting per-(branch,contingency) CSV rows so each row carries + // human-readable area/zone/owner names alongside the numbers. + std::map area_name_by_num; + std::map zone_name_by_num; + std::map owner_name_by_num; + { + boost::shared_ptr netdata = + pf_network->getNetworkData(); + int aT = 0, zT = 0, oT = 0; + netdata->getValue(AREA_TOTAL, &aT); + netdata->getValue(ZONE_TOTAL, &zT); + netdata->getValue(OWNER_TOTAL, &oT); + for (int i = 0; i < aT; i++) { + int n = 0; std::string s; + netdata->getValue(AREAINTG_NUMBER, &n, i); + netdata->getValue(AREAINTG_NAME, &s, i); + area_name_by_num[n] = s; + } + for (int i = 0; i < zT; i++) { + int n = 0; std::string s; + netdata->getValue(ZONE_NUMBER, &n, i); + netdata->getValue(ZONE_NAME, &s, i); + zone_name_by_num[n] = s; + } + for (int i = 0; i < oT; i++) { + int n = 0; std::string s; + netdata->getValue(OWNER_NUMBER, &n, i); + netdata->getValue(OWNER_NAME, &s, i); + owner_name_by_num[n] = s; + } + } + + // Per-rank bus metadata + wide/long branch-row outputs for csv_flat + // (long-form one row per branch per case) and csv_delta (wide-form one + // row per branch per case joining base and cont state). Both share the + // bus_meta load, the buses sidecar, and the per-rank .part-file gather. + bool wantBusSidecar = (outputFormat == "csv_flat" || + outputFormat == "csv_delta"); + struct BusMeta { + std::string name; + double basekv; + int area, zone, owner; + }; + std::map bus_meta; + if (wantBusSidecar) { + int nBus = pf_network->numBuses(); + for (int i = 0; i < nBus; i++) { + gridpack::powerflow::PFBus *bus = + dynamic_cast(pf_network->getBus(i).get()); + if (!bus) continue; + int orig = pf_network->getOriginalBusIndex(i); + BusMeta m; + m.name = bus->getBusName(); + m.basekv = bus->getBaseKV(); + m.area = bus->getArea(); + m.zone = bus->getZone(); + m.owner = bus->getOwner(); + bus_meta[orig] = m; + } + } + // (area, base kV) per bus, all-gathered over task_comm so filter lookups + // work on gathered rows. + struct BusAreaKv { int area; double basekv; }; + std::map bus_ak; + if (wantBusSidecar || haveMonitorFilter) { + std::vector lid, larea; + std::vector lkv; + int nBus = pf_network->numBuses(); + for (int i = 0; i < nBus; i++) { + if (!pf_network->getActiveBus(i)) continue; + gridpack::powerflow::PFBus *bus = + dynamic_cast(pf_network->getBus(i).get()); + if (!bus) continue; + lid.push_back(pf_network->getOriginalBusIndex(i)); + larea.push_back(bus->getArea()); + lkv.push_back(bus->getBaseKV()); + } + MPI_Comm tc = static_cast(task_comm); + int tsize = task_comm.size(); + int nloc = static_cast(lid.size()); + std::vector counts(tsize, 0), displs(tsize, 0); + MPI_Allgather(&nloc, 1, MPI_INT, &counts[0], 1, MPI_INT, tc); + int tot = 0; + for (int p = 0; p < tsize; p++) { displs[p] = tot; tot += counts[p]; } + std::vector gid(tot > 0 ? tot : 1), garea(tot > 0 ? tot : 1); + std::vector gkv(tot > 0 ? tot : 1); + int *lid_p = nloc > 0 ? &lid[0] : NULL; + int *larea_p = nloc > 0 ? &larea[0] : NULL; + double *lkv_p = nloc > 0 ? &lkv[0] : NULL; + MPI_Allgatherv(lid_p, nloc, MPI_INT, &gid[0], &counts[0], &displs[0], MPI_INT, tc); + MPI_Allgatherv(larea_p, nloc, MPI_INT, &garea[0], &counts[0], &displs[0], MPI_INT, tc); + MPI_Allgatherv(lkv_p, nloc, MPI_DOUBLE, &gkv[0], &counts[0], &displs[0], MPI_DOUBLE, tc); + for (int i = 0; i < tot; i++) { + BusAreaKv a; + a.area = garea[i]; + a.basekv = gkv[i]; + bus_ak[gid[i]] = a; + } + } + + // Strip surrounding single quotes (PSS/E style) and outer whitespace. + auto trim_quoted = [](const std::string &in) -> std::string { + std::string s = in; + size_t a = s.find_first_not_of(" \t"); + size_t b = s.find_last_not_of(" \t"); + if (a == std::string::npos) return std::string(); + s = s.substr(a, b - a + 1); + if (s.size() >= 2 && s.front() == '\'' && s.back() == '\'') { + s = s.substr(1, s.size() - 2); + } + a = s.find_first_not_of(" \t"); + b = s.find_last_not_of(" \t"); + return (a == std::string::npos) ? std::string() : s.substr(a, b - a + 1); + }; + auto lookup_name = [&](const std::map &m, int n) -> std::string { + std::map::const_iterator it = m.find(n); + if (it == m.end()) return std::string(); + return trim_quoted(it->second); + }; + + // Per-rank bus metadata sidecar (deduped by world rank 0 after the loop). + if (wantBusSidecar) { + std::ostringstream oss; + oss << outputFile << "_buses." << world.rank() << ".part"; + std::ofstream fbus(oss.str().c_str(), + std::ios::out | std::ios::trunc | std::ios::binary); + fbus << std::fixed; + for (std::map::const_iterator it = bus_meta.begin(); + it != bus_meta.end(); ++it) { + const BusMeta &m = it->second; + fbus << it->first << "," + << trim_quoted(m.name) << "," + << std::setprecision(2) << m.basekv << "," + << m.area << "," << m.zone << "," << m.owner << "," + << lookup_name(area_name_by_num, m.area) << "," + << lookup_name(zone_name_by_num, m.zone) << "," + << lookup_name(owner_name_by_num, m.owner) + << "\n"; + } + fbus.close(); + } + + // Per-rank streaming output file. Opened on first row written so non-csv_flat + // runs and ranks that produce no rows leave nothing behind. + std::string flatPartPath; + if (outputFormat == "csv_flat") { + std::ostringstream oss; + oss << outputFile << "_flat." << world.rank() << ".part"; + flatPartPath = oss.str(); + } + std::ofstream flatPart; + size_t flatRowCount = 0; + + // Per-rank _violations.csv stream. Populated by every output format that + // knows about violations (json/csv via populateViolations; csv_flat/csv_delta + // inline where loading_pct is already computed). Emit-branch and emit-voltage + // helpers below open the file lazily. + std::string violPartPath; + { + std::ostringstream oss; + oss << outputFile << "_violations." << world.rank() << ".part"; + violPartPath = oss.str(); + } + std::ofstream violPart; + size_t violRowCount = 0; + // Running summary state; captureFlatRows / captureDeltaRows / populateViolations + // all funnel through the emit helpers, so this stays consistent regardless of + // outputFormat. + struct WorstBranchState { + double loading_pct = 0.0; + int from = 0, to = 0; + std::string ckt; + std::string ct_name; + }; + struct WorstVoltageState { + double v_pu = 1.0; + double dev_pu = 0.0; // signed + int bus_id = 0; + std::string ct_name; + }; + WorstBranchState worstBr; + WorstVoltageState worstVLo; + worstVLo.v_pu = 1e9; // start high so any real low value beats it + WorstVoltageState worstVHi; + worstVHi.v_pu = -1e9; + std::set ctsWithBranchViol; + std::set ctsWithVoltageViol; + // Rank-local per-ct composite indices; reduced to world 0 in summary. + std::map ctPi; // sum(mva/rate)^2, all monitored branches + std::map ctVpi; // sum((v-1)/dv)^2, all energized buses + std::map ctVdev; // sum(v-limit)^2, violated buses only + std::map ctWorstLoading; // max loading_pct among violated branches + std::map ctWorstVdev; // max |v - limit| among violated buses + std::map ctWorstVpu; // v_pu that produced ctWorstVdev + auto accumBranchPi = [&](const std::string &ct_name, + double mva, double rate) { + if (rate <= 0.0) return; + double r = mva / rate; + ctPi[ct_name] += r * r; + }; + // Textbook voltage PI: normalize by the direction-aware half-band so a bus + // at Vmin or Vmax contributes 1.0 even when limits are asymmetric. + // Skip isolated (v<=0) and NaN/inf buses so numerical failures on a single + // bus don't poison the accumulator. + auto accumVoltagePi = [&](const std::string &ct_name, double v_pu) { + if (v_pu <= 0.0 || !std::isfinite(v_pu)) return; + double denom = (v_pu < 1.0) ? (1.0 - Vmin) : (Vmax - 1.0); + if (denom <= 0.0) return; + double r = (v_pu - 1.0) / denom; + ctVpi[ct_name] += r * r; + }; + // Legacy deviation metric: (v-limit)^2 accrued only for violated buses. + auto accumVoltageDev = [&](const std::string &ct_name, double dev_pu) { + ctVdev[ct_name] += dev_pu * dev_pu; + }; + auto openViolPart = [&]() { + if (violPart.is_open()) return; + violPart.open(violPartPath.c_str(), std::ios::out | std::ios::trunc); + violPart << std::fixed; + }; + // Row schema (same 11 columns for branch and voltage rows; unused fields blank): + // event_idx,contingency,type,element,mva_or_vpu,rate_or_limit,loading_percent, + // base_mva,delta,severity + // For branch: element = "from-to-ckt", mva_or_vpu = MVA, rate_or_limit = rate, + // loading_percent = %, base_mva/delta populated, severity. + // For voltage: element = "bus_id", mva_or_vpu = v_pu, rate_or_limit = "low/high", + // loading_percent = "", base_mva = "", delta = signed dev pu. + auto emitBranchViolation = [&](int event_idx, const std::string &ct_name, + int from, int to, const std::string &ckt, + double mva, double rate, + double loading_pct, double base_mva) { + openViolPart(); + const char *sev = (loading_pct >= 105.0) ? "critical" : "warning"; + violPart << event_idx << "," << ct_name << ",branch," + << from << "-" << to << "-" << ckt << "," + << std::setprecision(4) << mva << "," + << std::setprecision(4) << rate << "," + << std::setprecision(2) << loading_pct << "," + << std::setprecision(4) << base_mva << "," + << std::setprecision(4) << (mva - base_mva) << "," + << sev << "\n"; + violRowCount++; + if (loading_pct > worstBr.loading_pct) { + worstBr.loading_pct = loading_pct; + worstBr.from = from; worstBr.to = to; + worstBr.ckt = ckt; + worstBr.ct_name = ct_name; + } + ctsWithBranchViol.insert(ct_name); + double &pc = ctWorstLoading[ct_name]; + if (loading_pct > pc) pc = loading_pct; + }; + auto emitVoltageViolation = [&](int event_idx, const std::string &ct_name, + int bus_id, double v_pu, + double lo, double hi) { + openViolPart(); + double dev = (v_pu < lo) ? (v_pu - lo) + : (v_pu > hi) ? (v_pu - hi) + : 0.0; + if (dev == 0.0) return; + const char *sev = (std::abs(dev) >= 0.05) ? "critical" : "warning"; + // rate_or_limit column holds "low_pu:high_pu" for voltage rows. + std::ostringstream limits; + limits << std::setprecision(4) << std::fixed << lo << ":" << hi; + violPart << event_idx << "," << ct_name << ",voltage," + << bus_id << "," + << std::setprecision(6) << v_pu << "," + << limits.str() << "," + << "," // loading_percent blank + << "," // base_mva blank + << std::setprecision(6) << dev << "," + << sev << "\n"; + violRowCount++; + if (dev < 0.0 && v_pu < worstVLo.v_pu) { + worstVLo.v_pu = v_pu; + worstVLo.dev_pu = dev; + worstVLo.bus_id = bus_id; + worstVLo.ct_name = ct_name; + } + if (dev > 0.0 && v_pu > worstVHi.v_pu) { + worstVHi.v_pu = v_pu; + worstVHi.dev_pu = dev; + worstVHi.bus_id = bus_id; + worstVHi.ct_name = ct_name; + } + ctsWithVoltageViol.insert(ct_name); + accumVoltageDev(ct_name, dev); + double absDev = std::abs(dev); + double &d = ctWorstVdev[ct_name]; + if (absDev > d) { d = absDev; ctWorstVpu[ct_name] = v_pu; } + }; + + // (from, to, ckt) key shared by the monitor allowlist and base_cache. + struct BranchKey { + int from, to; + std::string ckt; + bool operator<(const BranchKey &o) const { + if (from != o.from) return from < o.from; + if (to != o.to ) return to < o.to; + return ckt < o.ckt; + } + }; + + // Per-branch rate-A/B/C from the parsed network data. Keyed by (from,to,ckt) + // so the csv_flat / csv_delta emit paths can pick the configured rating. + // Built once before the contingency loop. Each rank only sees its own + // active+ghost branches; that's fine -- the emit path is also rank-local. + struct BranchRates { + double rate_a, rate_b, rate_c; + }; + std::map branch_rates; + if (outputFormat == "csv_flat" || outputFormat == "csv_delta") { + int nBranch = pf_network->numBranches(); + for (int i = 0; i < nBranch; i++) { + boost::shared_ptr bd = + pf_network->getBranchData(i); + if (!bd) continue; + int from = 0, to = 0, nelems = 0; + bd->getValue(BRANCH_FROMBUS, &from); + bd->getValue(BRANCH_TOBUS, &to); + if (!bd->getValue(BRANCH_NUM_ELEMENTS, &nelems)) continue; + for (int k = 0; k < nelems; k++) { + std::string ckt; + if (!bd->getValue(BRANCH_CKT, &ckt, k)) continue; + // Trim leading/trailing whitespace and PSS/E surrounding quotes. + size_t a = ckt.find_first_not_of(" \t"); + size_t b = ckt.find_last_not_of(" \t"); + ckt = (a == std::string::npos) ? std::string() + : ckt.substr(a, b - a + 1); + if (ckt.size() >= 2 && ckt.front() == '\'' && ckt.back() == '\'') { + ckt = ckt.substr(1, ckt.size() - 2); + a = ckt.find_first_not_of(" \t"); + b = ckt.find_last_not_of(" \t"); + ckt = (a == std::string::npos) ? std::string() + : ckt.substr(a, b - a + 1); + } + BranchRates r; + r.rate_a = 0.0; r.rate_b = 0.0; r.rate_c = 0.0; + bd->getValue(BRANCH_RATING_A, &r.rate_a, k); + bd->getValue(BRANCH_RATING_B, &r.rate_b, k); + bd->getValue(BRANCH_RATING_C, &r.rate_c, k); + BranchKey key; + key.from = from; key.to = to; key.ckt = ckt; + branch_rates[key] = r; + } + } + } + // Base case always uses rate-A (PSS/E "normal" rating). Contingency rows + // use whichever the user picked, with A->B->C fallback if zero/missing. + auto pickContRate = [&](const BranchRates &r) -> double { + if (contingencyRating == "A") { + return r.rate_a; + } + if (contingencyRating == "B") { + return (r.rate_b > 0.0) ? r.rate_b : r.rate_a; + } + if (r.rate_c > 0.0) return r.rate_c; + if (r.rate_b > 0.0) return r.rate_b; + return r.rate_a; + }; + + // Monitor allowlist parsed from monitorBranchesFile. Empty -> emit all. + std::set monitorSet; + if (!monitorBranchesFile.empty()) { + std::ifstream fin(monitorBranchesFile.c_str()); + if (!fin.is_open()) { + if (world.rank() == 0) { + printf("WARNING: monitorBranchesFile '%s' not found; emitting all branches\n", + monitorBranchesFile.c_str()); + } + } else { + std::string line; + size_t lineNo = 0; + while (std::getline(fin, line)) { + lineNo++; + // Strip trailing CR (Windows line endings). + while (!line.empty() && (line[line.size()-1] == '\r' || + line[line.size()-1] == '\n')) { + line.resize(line.size()-1); + } + // Skip blank lines and comments. + size_t firstNon = line.find_first_not_of(" \t"); + if (firstNon == std::string::npos) continue; + if (line[firstNon] == '#') continue; + // Tokenize on commas. + std::vector tok; + size_t pos = 0; + while (pos <= line.size()) { + size_t comma = line.find(',', pos); + std::string t = (comma == std::string::npos) + ? line.substr(pos) + : line.substr(pos, comma - pos); + size_t a = t.find_first_not_of(" \t"); + size_t b = t.find_last_not_of(" \t"); + tok.push_back((a == std::string::npos) ? std::string() + : t.substr(a, b - a + 1)); + if (comma == std::string::npos) break; + pos = comma + 1; + } + if (tok.size() < 3) continue; + // Skip header row: any non-numeric first field. + if (tok[0].empty()) continue; + bool numeric = true; + for (size_t ci = 0; ci < tok[0].size(); ci++) { + char c = tok[0][ci]; + if (!(c >= '0' && c <= '9') && c != '-' && c != '+') { + numeric = false; break; + } + } + if (!numeric) continue; + BranchKey k; + k.from = atoi(tok[0].c_str()); + k.to = atoi(tok[1].c_str()); + k.ckt = tok[2]; + monitorSet.insert(k); + } + if (world.rank() == 0) { + printf("Monitor allowlist: %zu branches loaded from %s\n", + monitorSet.size(), monitorBranchesFile.c_str()); + } + } + } + // Endpoints of allowlisted branches; gates voltage rows under an allowlist. + std::set monitorBusSet; + for (std::set::const_iterator it = monitorSet.begin(); + it != monitorSet.end(); ++it) { + monitorBusSet.insert(it->from); + monitorBusSet.insert(it->to); + } + // Area/kV gate. Either-endpoint match for areas (catches tie-lines). + // kV is gated on max(kv_from, kv_to) so a 138/13.8 stepdown counts as 138. + // Empty area set / zero kV bound = unrestricted on that dimension. + auto passesAreaKv = [&](int area_from, int area_to, + double kv_from, double kv_to) { + if (!monitorAreas.empty()) { + if (monitorAreas.find(area_from) == monitorAreas.end() && + monitorAreas.find(area_to) == monitorAreas.end()) { + return false; + } + } + double kv_max = (kv_from > kv_to) ? kv_from : kv_to; + if (monitorKvMin > 0.0 && kv_max < monitorKvMin) return false; + if (monitorKvMax > 0.0 && kv_max > monitorKvMax) return false; + return true; + }; + // When monitorBranchesFile presents, it overrides area/kV criteria. + bool haveAreaKvFilter = !monitorAreas.empty() || + monitorKvMin > 0.0 || + monitorKvMax > 0.0; + if (!monitorSet.empty() && haveAreaKvFilter) { + if (world.rank() == 0) { + printf("WARNING: monitorBranchesFile is set; ignoring " + "monitorAreas/monitorKvMin/monitorKvMax\n"); + } + monitorAreas.clear(); + monitorKvMin = 0.0; + monitorKvMax = 0.0; + haveAreaKvFilter = false; + } + if (world.rank() == 0) { + if (!monitorAreas.empty()) { + printf("Monitor areas filter: %zu areas\n", monitorAreas.size()); + } + if (monitorKvMin > 0.0 || monitorKvMax > 0.0) { + printf("Monitor kV filter: min=%.2f max=%.2f (0 means unbounded)\n", + monitorKvMin, monitorKvMax); + } + } + auto busAreaKv = [&](int bus, int &area, double &kv) { + std::map::const_iterator it = bus_ak.find(bus); + area = (it != bus_ak.end()) ? it->second.area : 0; + kv = (it != bus_ak.end()) ? it->second.basekv : 0.0; + }; + // Branch monitor predicate shared by every output path (ckt padding trimmed). + auto branchMonitored = [&](int from, int to, const std::string &ckt) -> bool { + if (!monitorSet.empty()) { + BranchKey k; + k.from = from; k.to = to; k.ckt = ckt; + while (!k.ckt.empty() && (k.ckt[k.ckt.size()-1] == ' ' || + k.ckt[k.ckt.size()-1] == '\t')) + k.ckt.resize(k.ckt.size()-1); + return monitorSet.find(k) != monitorSet.end(); + } + if (!haveAreaKvFilter) return true; + int af = 0, at = 0; + double kf = 0.0, kt = 0.0; + busAreaKv(from, af, kf); + busAreaKv(to, at, kt); + return passesAreaKv(af, at, kf, kt); + }; + // Bus counterpart: allowlist endpoint, else own area / base kV. + auto busMonitored = [&](int bus) -> bool { + if (!monitorSet.empty()) { + return monitorBusSet.find(bus) != monitorBusSet.end(); + } + if (!haveAreaKvFilter) return true; + int area = 0; + double kv = 0.0; + busAreaKv(bus, area, kv); + if (!monitorAreas.empty() && monitorAreas.find(area) == monitorAreas.end()) + return false; + if (monitorKvMin > 0.0 && kv < monitorKvMin) return false; + if (monitorKvMax > 0.0 && kv > monitorKvMax) return false; + return true; + }; + + struct BaseFlow { + double p_mw, q_mvar, mva, loading_pct; + double base_rate, cont_rate; + double v_from_pu, v_to_pu, ang_from_deg, ang_to_deg; + double base_kv_from, base_kv_to; + int area_from, area_to; + }; + std::map base_cache; + std::string deltaPartPath; + if (outputFormat == "csv_delta") { + std::ostringstream oss; + oss << outputFile << "_delta." << world.rank() << ".part"; + deltaPartPath = oss.str(); + } + std::ofstream deltaPart; + size_t deltaRowCount = 0; + size_t deltaSkipCount = 0; + + // Convergence sidecar rows. + struct ConvRow { + int event_idx; + std::string name; + std::string type; + gridpack::utility::ConvergenceSummary cs; + std::string status; + }; + std::vector localConvRows; + // _convergence.csv: written for every outputFormat. + bool emitConv = true; + + // Lambda: parse current solved flow_str/vr_str and stream one CSV row + // per branch into the rank's .part file. Called once per converged case + // (base + each contingency) on every task communicator; non-rank-0 + // task_comm members short-circuit after the collective. + auto captureFlatRows = [&](int event_idx, const std::string &name, + bool emit, bool is_base) { + std::vector v_strs = pf_app.writeBusString("vr_str"); + std::vector b_strs = pf_app.writeBranchString("flow_str"); + if (!emit || task_comm.rank() != 0) return; + if (!flatPart.is_open()) { + flatPart.open(flatPartPath.c_str(), std::ios::out | std::ios::trunc); + flatPart << std::fixed; + } + std::map > vbymag_ang; + for (size_t vi = 0; vi < v_strs.size(); vi++) { + int bus_id = 0, use_vmag = 0, changed = 0; + double angle = 0.0, vmag = 0.0; + if (sscanf(v_strs[vi].c_str(), "%d %lf %lf %d %d", + &bus_id, &angle, &vmag, &use_vmag, &changed) == 5) { + vbymag_ang[bus_id] = std::make_pair(vmag, angle); + } + } + char ct_name[24]; + std::strncpy(ct_name, name.c_str(), sizeof(ct_name) - 1); + ct_name[sizeof(ct_name) - 1] = '\0'; + // Voltage PI and violations on monitored buses, contingency rows only. + if (!is_base) { + for (std::map >::const_iterator vit = + vbymag_ang.begin(); vit != vbymag_ang.end(); ++vit) { + if (!busMonitored(vit->first)) continue; + double v_pu = vit->second.first; + accumVoltagePi(ct_name, v_pu); + if (v_pu <= 0.0 || !std::isfinite(v_pu)) continue; + if (v_pu < Vmin || v_pu > Vmax) { + emitVoltageViolation(event_idx, ct_name, vit->first, v_pu, Vmin, Vmax); + } + } + } + for (size_t bi = 0; bi < b_strs.size(); bi++) { + char ckt_buf[16] = {0}; + int viol = 0; + double p = 0.0, q = 0.0, perf = 0.0, ratea = 0.0; + int from = 0, to = 0; + if (sscanf(b_strs[bi].c_str(), + "%d %d %15s %lf %lf %lf %lf %d", + &from, &to, ckt_buf, &p, &q, &perf, &ratea, &viol) != 8) { + continue; + } + char ckt[4]; + std::strncpy(ckt, ckt_buf, 3); ckt[3] = '\0'; + BranchKey mk; + mk.from = from; mk.to = to; mk.ckt = ckt; + while (!mk.ckt.empty() && mk.ckt[mk.ckt.size()-1] == ' ') + mk.ckt.resize(mk.ckt.size()-1); + if (!branchMonitored(from, to, mk.ckt)) continue; + std::map::const_iterator rIt = branch_rates.find(mk); + double rate_sel = ratea; + if (rIt != branch_rates.end()) { + rate_sel = is_base ? rIt->second.rate_a : pickContRate(rIt->second); + } + double flow_mva = std::sqrt(p*p + q*q); + double loading_pct = (rate_sel > 0.0) ? (flow_mva / rate_sel) * 100.0 : 0.0; + if (!is_base) accumBranchPi(ct_name, flow_mva, rate_sel); + // Stream to _violations.csv for csv_flat runs (contingency rows only). + if (!is_base && loading_pct > violationSeverityThreshold * 100.0) { + double base_mva = 0.0; + // No base_cache in csv_flat mode; look up in the persistent map built + // by populateBaseCache-style capture below? We don't have one for + // csv_flat, so base_mva stays 0 -- delta will just equal mva. + emitBranchViolation(event_idx, ct_name, from, to, ckt, + flow_mva, rate_sel, loading_pct, base_mva); + } + std::map >::const_iterator vf = + vbymag_ang.find(from); + std::map >::const_iterator vt = + vbymag_ang.find(to); + double v_from = (vf != vbymag_ang.end()) ? vf->second.first : 0.0; + double ang_from_deg = (vf != vbymag_ang.end()) ? vf->second.second : 0.0; + double v_to = (vt != vbymag_ang.end()) ? vt->second.first : 0.0; + double ang_to_deg = (vt != vbymag_ang.end()) ? vt->second.second : 0.0; + flatPart << event_idx << "," << ct_name << "," + << from << "," << to << "," << ckt << "," + << std::setprecision(4) << p << "," + << std::setprecision(4) << q << "," + << std::setprecision(4) << flow_mva << "," + << std::setprecision(4) << rate_sel << "," + << std::setprecision(2) << loading_pct << "," + << viol << "," + << std::setprecision(6) << v_from << "," + << std::setprecision(6) << v_to << "," + << std::setprecision(4) << ang_from_deg << "," + << std::setprecision(4) << ang_to_deg + << "\n"; + flatRowCount++; + } + }; + + // Populate base_cache from current solved state. Called once after base + // solve on every rank (csv_delta only); world.rank() == 0 is not special + // here -- each rank caches the branches it sees on its task_comm so it + // can join later in captureDeltaRows. + auto populateBaseCache = [&]() { + std::vector v_strs = pf_app.writeBusString("vr_str"); + std::vector b_strs = pf_app.writeBranchString("flow_str"); + if (task_comm.rank() != 0) return; + std::map > vbymag_ang; + for (size_t vi = 0; vi < v_strs.size(); vi++) { + int bus_id = 0, use_vmag = 0, changed = 0; + double angle = 0.0, vmag = 0.0; + if (sscanf(v_strs[vi].c_str(), "%d %lf %lf %d %d", + &bus_id, &angle, &vmag, &use_vmag, &changed) == 5) { + vbymag_ang[bus_id] = std::make_pair(vmag, angle); + } + } + for (size_t bi = 0; bi < b_strs.size(); bi++) { + char ckt_buf[16] = {0}; + int viol = 0; + double p = 0.0, q = 0.0, perf = 0.0, ratea = 0.0; + int from = 0, to = 0; + if (sscanf(b_strs[bi].c_str(), + "%d %d %15s %lf %lf %lf %lf %d", + &from, &to, ckt_buf, &p, &q, &perf, &ratea, &viol) != 8) { + continue; + } + BranchKey k; + k.from = from; k.to = to; + k.ckt = std::string(ckt_buf); + // Strip trailing spaces from ckt so the key matches what flow_str + // returns later (sscanf %15s already trims leading whitespace). + while (!k.ckt.empty() && k.ckt[k.ckt.size()-1] == ' ') k.ckt.resize(k.ckt.size()-1); + if (!branchMonitored(from, to, k.ckt)) continue; + double base_rate = ratea, cont_rate = ratea; + std::map::const_iterator rIt = branch_rates.find(k); + if (rIt != branch_rates.end()) { + base_rate = rIt->second.rate_a; + cont_rate = pickContRate(rIt->second); + } + BaseFlow bf; + bf.p_mw = p; + bf.q_mvar = q; + bf.mva = std::sqrt(p*p + q*q); + bf.base_rate = base_rate; + bf.cont_rate = cont_rate; + bf.loading_pct = (base_rate > 0.0) ? (bf.mva / base_rate) * 100.0 : 0.0; + std::map >::const_iterator vf = + vbymag_ang.find(from); + std::map >::const_iterator vt = + vbymag_ang.find(to); + bf.v_from_pu = (vf != vbymag_ang.end()) ? vf->second.first : 0.0; + bf.ang_from_deg = (vf != vbymag_ang.end()) ? vf->second.second : 0.0; + bf.v_to_pu = (vt != vbymag_ang.end()) ? vt->second.first : 0.0; + bf.ang_to_deg = (vt != vbymag_ang.end()) ? vt->second.second : 0.0; + busAreaKv(from, bf.area_from, bf.base_kv_from); + busAreaKv(to, bf.area_to, bf.base_kv_to); + base_cache[k] = bf; + } + }; + + // Wide-form (base+cont on same row) capture for csv_delta. Mirrors + // captureFlatRows but joins each branch with base_cache. Branches not + // in base_cache are counted in deltaSkipCount and skipped silently. + auto captureDeltaRows = [&](int event_idx, + const gridpack::powerflow::Contingency &evt, + bool emit) { + std::vector v_strs = pf_app.writeBusString("vr_str"); + std::vector b_strs = pf_app.writeBranchString("flow_str"); + if (!emit || task_comm.rank() != 0) return; + if (!deltaPart.is_open()) { + deltaPart.open(deltaPartPath.c_str(), std::ios::out | std::ios::trunc); + deltaPart << std::fixed; + } + // cont_event_facility: disabled -- join event_idx against + // _contingencies.csv instead, which names every outaged element + // rather than the first plus "(+N more)". Kept commented in case the column + // is wanted back; uncomment this block, the emission below and the header + // field together. The [area] lookup needs a complete bus_meta (groupSize=1). + // std::string facility; + // // clean2Char pads ids to two chars; trim so the label has no stray space. + // auto rtrimId = [](const std::string &in) -> std::string { + // std::string t = in; + // while (!t.empty() && (t[t.size()-1] == ' ' || t[t.size()-1] == '\t')) + // t.resize(t.size()-1); + // return t; + // }; + // // "[area] " for both kinds; the type column says which. + // if (evt.p_type == Branch && !evt.p_from.empty()) { + // int outFrom = evt.p_from[0]; + // int area = 0; + // std::map::const_iterator mf = bus_meta.find(outFrom); + // if (mf != bus_meta.end()) area = mf->second.area; + // char buf[64]; + // snprintf(buf, sizeof(buf), "[%d] %d %d %s", + // area, outFrom, evt.p_to[0], rtrimId(evt.p_ckt[0]).c_str()); + // facility = buf; + // if (evt.p_from.size() > 1) { + // char suf[24]; + // snprintf(suf, sizeof(suf), " (+%zu more)", evt.p_from.size() - 1); + // facility += suf; + // } + // } else if (evt.p_type == Generator && !evt.p_busid.empty()) { + // int outBus = evt.p_busid[0]; + // int area = 0; + // std::map::const_iterator mg = bus_meta.find(outBus); + // if (mg != bus_meta.end()) area = mg->second.area; + // char buf[64]; + // snprintf(buf, sizeof(buf), "[%d] %d %s", + // area, outBus, rtrimId(evt.p_genid[0]).c_str()); + // facility = buf; + // if (evt.p_busid.size() > 1) { + // char suf[24]; + // snprintf(suf, sizeof(suf), " (+%zu more)", evt.p_busid.size() - 1); + // facility += suf; + // } + // } + std::string ct_name = evt.p_name; + while (!ct_name.empty() && ct_name[ct_name.size()-1] == ' ') + ct_name.resize(ct_name.size()-1); + const char *type_str = (evt.p_type == Branch) ? "branch" : "generator"; + std::map > vbymag_ang; + for (size_t vi = 0; vi < v_strs.size(); vi++) { + int bus_id = 0, use_vmag = 0, changed = 0; + double angle = 0.0, vmag = 0.0; + if (sscanf(v_strs[vi].c_str(), "%d %lf %lf %d %d", + &bus_id, &angle, &vmag, &use_vmag, &changed) == 5) { + vbymag_ang[bus_id] = std::make_pair(vmag, angle); + } + } + // Voltage PI and violations on monitored buses. + for (std::map >::const_iterator vit = + vbymag_ang.begin(); vit != vbymag_ang.end(); ++vit) { + if (!busMonitored(vit->first)) continue; + double v_pu = vit->second.first; + accumVoltagePi(ct_name, v_pu); + if (v_pu <= 0.0 || !std::isfinite(v_pu)) continue; + if (v_pu < Vmin || v_pu > Vmax) { + emitVoltageViolation(event_idx, ct_name, vit->first, v_pu, Vmin, Vmax); + } + } + for (size_t bi = 0; bi < b_strs.size(); bi++) { + char ckt_buf[16] = {0}; + int viol = 0; + double p = 0.0, q = 0.0, perf = 0.0, ratea = 0.0; + int from = 0, to = 0; + if (sscanf(b_strs[bi].c_str(), + "%d %d %15s %lf %lf %lf %lf %d", + &from, &to, ckt_buf, &p, &q, &perf, &ratea, &viol) != 8) { + continue; + } + BranchKey k; + k.from = from; k.to = to; + k.ckt = std::string(ckt_buf); + while (!k.ckt.empty() && k.ckt[k.ckt.size()-1] == ' ') + k.ckt.resize(k.ckt.size()-1); + if (!branchMonitored(from, to, k.ckt)) continue; + std::map::const_iterator it = base_cache.find(k); + if (it == base_cache.end()) { deltaSkipCount++; continue; } + const BaseFlow &bf = it->second; + double cont_mva = std::sqrt(p*p + q*q); + double cont_loading = (bf.cont_rate > 0.0) ? (cont_mva / bf.cont_rate) * 100.0 : 0.0; + accumBranchPi(ct_name, cont_mva, bf.cont_rate); + // Stream to _violations.csv (delta path knows base_mva already). + if (cont_loading > violationSeverityThreshold * 100.0) { + emitBranchViolation(event_idx, ct_name, from, to, k.ckt, + cont_mva, bf.cont_rate, cont_loading, bf.mva); + } + std::map >::const_iterator vf = + vbymag_ang.find(from); + std::map >::const_iterator vt = + vbymag_ang.find(to); + double v_from_c = (vf != vbymag_ang.end()) ? vf->second.first : 0.0; + double a_from_c = (vf != vbymag_ang.end()) ? vf->second.second : 0.0; + double v_to_c = (vt != vbymag_ang.end()) ? vt->second.first : 0.0; + double a_to_c = (vt != vbymag_ang.end()) ? vt->second.second : 0.0; + double d_ang_b = bf.ang_from_deg - bf.ang_to_deg; + double d_ang_c = a_from_c - a_to_c; + // Across-branch drop, same convention as the angle deltas above. + double d_v_b = bf.v_from_pu - bf.v_to_pu; + double d_v_c = v_from_c - v_to_c; + deltaPart << event_idx << "," << ct_name << "," << type_str << "," + << from << "," << to << "," << k.ckt << "," + << std::setprecision(2) << bf.base_kv_from << "," + << std::setprecision(2) << bf.base_kv_to << "," + << bf.area_from << "," << bf.area_to << "," + << std::setprecision(4) << bf.base_rate << "," + << std::setprecision(4) << bf.cont_rate << "," + << std::setprecision(4) << bf.p_mw << "," + << std::setprecision(4) << p << "," + << std::setprecision(4) << bf.q_mvar << "," + << std::setprecision(4) << q << "," + << std::setprecision(4) << bf.mva << "," + << std::setprecision(4) << cont_mva << "," + << std::setprecision(2) << bf.loading_pct << "," + << std::setprecision(2) << cont_loading << "," + << std::setprecision(6) << bf.v_from_pu << "," + << std::setprecision(6) << v_from_c << "," + << std::setprecision(6) << bf.v_to_pu << "," + << std::setprecision(6) << v_to_c << "," + << std::setprecision(4) << bf.ang_from_deg << "," + << std::setprecision(4) << a_from_c << "," + << std::setprecision(4) << bf.ang_to_deg << "," + << std::setprecision(4) << a_to_c << "," + << std::setprecision(6) << d_v_b << "," + << std::setprecision(6) << d_v_c << "," + << std::setprecision(4) << d_ang_b << "," + << std::setprecision(4) << d_ang_c + // << "," << facility // cont_event_facility: disabled + << "\n"; + deltaRowCount++; + } + }; + // Set minimum and maximum voltage limits on all buses pf_app.setVoltageLimits(Vmin, Vmax); - // Solve the base power flow calculation. This calculation is replicated on - // all task communicators - pf_app.solve(); - // Check for Qlimit violations - if (check_Qlim && !pf_app.checkQlimViolations()) { - pf_app.solve(); + // Route CA violation checks and loadingPercent through the same rating tier. + pf_app.setContingencyRating(contingencyRating); + // Solve the base power flow on every task communicator. Abort if it fails. + bool baseSolveOk = false; + try { + baseSolveOk = pf_app.solve(); + if (baseSolveOk && check_Qlim && !pf_app.checkQlimViolations()) { + baseSolveOk = pf_app.solve(); + } + } catch (const std::exception &e) { + if (world.rank() == 0) { + printf("ERROR: base-case solve threw exception: %s\n", e.what()); + } + baseSolveOk = false; + } catch (...) { + if (world.rank() == 0) { + printf("ERROR: base-case solve threw unknown exception\n"); + } + baseSolveOk = false; } - // Some buses may violate the voltage limits in the base problem. Flag these - // buses to ignore voltage violations on them. + if (!baseSolveOk) { + if (world.rank() == 0) { + gridpack::utility::ConvergenceSummary cs = pf_app.getConvergence(); + printf("ERROR: base case did not converge " + "(iterations=%d, final_tol=%.6e, " + "max_p_bus=%d max_p_mismatch=%.4f, " + "max_q_bus=%d max_q_mismatch=%.4f). " + "Aborting contingency analysis.\n", + cs.iterations, cs.finalTolerance, + cs.finalMismatch.maxPBus, cs.finalMismatch.maxPMismatch, + cs.finalMismatch.maxQBus, cs.finalMismatch.maxQMismatch); + } + world.barrier(); + MPI_Abort(static_cast(world), 1); + } + // Suppress voltage violations already present at base. pf_app.ignoreVoltageViolations(); - // Collect base case results for export + // Flag non-monitored elements "ignore" so the violation checks and .out + // listings follow the filter. + if (haveMonitorFilter) { + long cnt[4] = { 0, 0, 0, 0 }; // mon buses, buses, mon elems, elems + int nBus = pf_network->numBuses(); + for (int i = 0; i < nBus; i++) { + gridpack::powerflow::PFBus *bus = + dynamic_cast(pf_network->getBus(i).get()); + if (!bus) continue; + bool mon = busMonitored(pf_network->getOriginalBusIndex(i)); + if (!mon) bus->setIgnore(true); + if (pf_network->getActiveBus(i)) { cnt[1]++; if (mon) cnt[0]++; } + } + int nBranch = pf_network->numBranches(); + for (int i = 0; i < nBranch; i++) { + gridpack::powerflow::PFBranch *br = + dynamic_cast(pf_network->getBranch(i).get()); + if (!br) continue; + int from = br->getBus1OriginalIndex(); + int to = br->getBus2OriginalIndex(); + std::vector tags = br->getLineTags(); + bool active = pf_network->getActiveBranch(i); + for (size_t t = 0; t < tags.size(); t++) { + bool mon = branchMonitored(from, to, tags[t]); + if (!mon) br->setIgnore(tags[t], true); + if (active) { cnt[3]++; if (mon) cnt[2]++; } + } + } + long tot[4] = { 0, 0, 0, 0 }; + MPI_Allreduce(cnt, tot, 4, MPI_LONG, MPI_SUM, + static_cast(task_comm)); + if (world.rank() == 0) { + printf("Monitor filter: %ld of %ld branch elements and %ld of %ld buses " + "monitored\n", tot[2], tot[3], tot[0], tot[1]); + if (tot[2] == 0 && tot[0] == 0) { + printf("WARNING: monitor filter matches nothing; all outputs will be " + "empty. Check monitorAreas/monitorKvMin/monitorKvMax/" + "monitorBranchesFile against the case.\n"); + } + } + } + // Drop non-monitored elements from a collected result set. + auto filterResults = [&](gridpack::utility::PowerFlowResults &r) { + if (!haveMonitorFilter) return; + std::vector buses; + for (size_t i = 0; i < r.buses.size(); i++) { + if (busMonitored(r.buses[i].busId)) buses.push_back(r.buses[i]); + } + r.buses.swap(buses); + std::vector branches; + for (size_t i = 0; i < r.branches.size(); i++) { + const gridpack::utility::BranchResult &b = r.branches[i]; + if (branchMonitored(b.fromBus, b.toBus, b.circuitId)) branches.push_back(b); + } + r.branches.swap(branches); + std::vector gens; + for (size_t i = 0; i < r.generators.size(); i++) { + if (busMonitored(r.generators[i].busId)) gens.push_back(r.generators[i]); + } + r.generators.swap(gens); + }; + + // Collect base case results for export. csv_flat captures rows directly + // in the hot loop and skips the heavyweight collectResults() path. gridpack::utility::PowerFlowResults baseCaseResults; - if (outputFormat != "text") { + if (outputFormat == "json" || outputFormat == "csv" || + outputFormat == "text") { baseCaseResults = pf_app.collectResults(); + filterResults(baseCaseResults); + } + if (outputFormat == "csv_flat") { + // The base case is replicated on every task communicator. captureFlatRows + // calls writeBusString/writeBranchString which are task_comm collectives, + // so every task_comm participates -- but only world rank 0 emits rows so + // the base case isn't duplicated in the final file. + captureFlatRows(0, std::string("base_case"), world.rank() == 0, true); + } + if (outputFormat == "csv_delta") { + // Cache base-case branch state on every rank for the contingency join. + populateBaseCache(); } // Check if auto-generation of N-1 contingencies is enabled @@ -545,8 +1595,55 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) printf("==================================================================\n\n"); } - // Print contingency details + // event_idx lookup table: one row per contingency (0 = base case); N-k + // element ids are ';'-separated within the columns. if (world.rank() == 0) { + std::string ctgFile = outputFile + "_contingencies.csv"; + std::ofstream cout_ctg(ctgFile.c_str(), std::ios::out | std::ios::trunc); + cout_ctg << "event_idx,contingency,type,n_elements," + "from_bus,to_bus,circuit_id,gen_bus,gen_id\n"; + size_t ctgRows = 0; + cout_ctg << "0,base_case,base,0,,,,,\n"; + ctgRows++; + // Trim clean2Char padding so ids join against the other CSVs. + auto rtrim = [](const std::string &in) -> std::string { + std::string t = in; + while (!t.empty() && (t[t.size()-1] == ' ' || t[t.size()-1] == '\t')) + t.resize(t.size()-1); + return t; + }; + for (size_t ei = 0; ei < events.size(); ei++) { + const gridpack::powerflow::Contingency &e = events[ei]; + int event_idx = static_cast(ei) + 1; + std::string nm = rtrim(e.p_name); + size_t n = 0; + const char *ty = "unknown"; + if (e.p_type == Branch) { n = e.p_from.size(); ty = "branch"; } + else if (e.p_type == Generator) { n = e.p_busid.size(); ty = "generator"; } + std::ostringstream c_from, c_to, c_ckt, c_gbus, c_gid; + for (size_t j = 0; j < n; j++) { + const char *sep = (j > 0) ? ";" : ""; + if (e.p_type == Branch) { + c_from << sep << e.p_from[j]; + c_to << sep << e.p_to[j]; + c_ckt << sep << rtrim(e.p_ckt[j]); + } else { + c_gbus << sep << e.p_busid[j]; + c_gid << sep << rtrim(e.p_genid[j]); + } + } + // Empty events still get a row so every event_idx decodes. + cout_ctg << event_idx << "," << nm << "," << ty << "," << n << "," + << c_from.str() << "," << c_to.str() << "," << c_ckt.str() << "," + << c_gbus.str() << "," << c_gid.str() << "\n"; + ctgRows++; + } + cout_ctg.close(); + printf("[contingencies] wrote %zu rows to %s\n", ctgRows, ctgFile.c_str()); + } + + // Print contingency details (gated on printCalcFiles; noisy for large lists) + if (print_calcs && world.rank() == 0) { int idx; for (idx = 0; idx < events.size(); idx++) { printf("Name: %s\n",events[idx].p_name.c_str()); @@ -579,164 +1676,158 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) int nbus = pf_network->totalBuses(); // Get bus voltage information for base case int i, j; -#ifdef USE_STATBLOCK + // StatBlock objects and the per-case scratch vectors live across the + // contingency loop so they are declared up here, regardless of whether + // statistics output is enabled. + boost::scoped_ptr vmag_stats; + boost::scoped_ptr vang_stats; + boost::scoped_ptr pgen_stats; + boost::scoped_ptr qgen_stats; + boost::scoped_ptr pflow_stats; + boost::scoped_ptr qflow_stats; + boost::scoped_ptr perf_stats; + std::vector v_vals; + int nsize = 0; + std::vector vmag, vang, pgen, qgen, pflow, qflow, perf; + std::vector mask, mag_mask; int t_store = timer->createCategory("Store Statistics"); - timer->start(t_store); - std::vector v_vals = pf_app.writeBusString("vr_str"); - int nsize = v_vals.size(); - std::vector mag_ids; - std::vector ids; - std::vector branch_ids; - std::vector mag_tags; - std::vector tags; - std::vector vmag; - std::vector vang; - std::vector mag_mask; - std::vector mask; - // Find bus IDs and create a dummy tag label and get voltage magnitude - // and angle for base case - for (i=0; i tokens = util.blankTokenizer(v_vals[i]); - int not_isolated = atoi(tokens[3].c_str()); - if (not_isolated == 1) { - mag_ids.push_back(atoi(tokens[0].c_str())); - mag_tags.push_back("1 "); - vmag.push_back(atof(tokens[2].c_str())); - if (atoi(tokens[4].c_str()) != 0) { - mag_mask.push_back(2); - } else { - mag_mask.push_back(1); + if (write_stats) { + timer->start(t_store); + v_vals = pf_app.writeBusString("vr_str"); + nsize = v_vals.size(); + std::vector mag_ids; + std::vector ids; + std::vector mag_tags; + std::vector tags; + // Find bus IDs and create a dummy tag label and get voltage magnitude + // and angle for base case + for (i=0; i tokens = util.blankTokenizer(v_vals[i]); + if (!busMonitored(atoi(tokens[0].c_str()))) continue; + int not_isolated = atoi(tokens[3].c_str()); + if (not_isolated == 1) { + mag_ids.push_back(atoi(tokens[0].c_str())); + mag_tags.push_back("1 "); + vmag.push_back(atof(tokens[2].c_str())); + if (atoi(tokens[4].c_str()) != 0) { + mag_mask.push_back(2); + } else { + mag_mask.push_back(1); + } } - } - ids.push_back(atoi(tokens[0].c_str())); - tags.push_back("1 "); - vang.push_back(atof(tokens[1].c_str())); - mask.push_back(1); - } - int nmags = vmag.size(); - world.max(&nmags,1); - world.max(&nbus,1); -#endif - // Create StatBlock objects for voltage magnitude and angles and add - // bus IDs to it -#ifdef USE_STATBLOCK - gridpack::analysis::StatBlock vmag_stats(world,nmags,ntasks+1); - gridpack::analysis::StatBlock vang_stats(world,nbus,ntasks+1); -#endif - // Add bus IDs and tags to StatBlock objects as well as base case values of - // voltage magnitude and angle -#ifdef USE_STATBLOCK - if (world.rank() == 0) { - vmag_stats.addRowLabels(mag_ids, mag_tags); - vang_stats.addRowLabels(ids, tags); - vmag_stats.addColumnValues(0,vmag,mag_mask); - vang_stats.addColumnValues(0,vang,mask); - } -#endif - // Get generator power information -#ifdef USE_STATBLOCK - v_vals.clear(); - ids.clear(); - tags.clear(); - mask.clear(); - std::vector pgen; - std::vector qgen; - v_vals = pf_app.writeBusString("power"); - nsize = v_vals.size(); - // Find bus IDs and tags for generators and eveluate Pg and Qg for base case - for (i=0; i tokens = util.blankTokenizer(v_vals[i]); - if (tokens.size()%4 != 0) { - printf("Incorrect generator listing\n"); - continue; - } - int ngen = tokens.size()/4; - for (j=0; j id1; - std::vector id2; - std::vector pmin, pmax; - std::vector pflow; - std::vector qflow; - std::vector perf; - v_vals = pf_app.writeBranchString("flow_str"); - nsize = v_vals.size(); - // Parse branch line endpoints as well as line IDs and values of P and Q for - // base case - for (i=0; i tokens = util.blankTokenizer(v_vals[i]); - if (tokens.size()%8 != 0) { - printf("Incorrect branch power flow listing\n"); - continue; + int nmags = vmag.size(); + world.max(&nmags,1); + world.max(&nbus,1); + // Create StatBlock objects for voltage magnitude and angles and add + // bus IDs to it as well as base case values + vmag_stats.reset(new gridpack::analysis::StatBlock(world,nmags,ntasks+1)); + vang_stats.reset(new gridpack::analysis::StatBlock(world,nbus,ntasks+1)); + if (world.rank() == 0) { + vmag_stats->addRowLabels(mag_ids, mag_tags); + vang_stats->addRowLabels(ids, tags); + vmag_stats->addColumnValues(0,vmag,mag_mask); + vang_stats->addColumnValues(0,vang,mask); } - int nline = tokens.size()/8; - for (j=0; j tokens = util.blankTokenizer(v_vals[i]); + if (tokens.size()%4 != 0) { + printf("Incorrect generator listing\n"); + continue; + } + int ngen = tokens.size()/4; + for (j=0; jaddRowLabels(ids, tags); + qgen_stats->addRowLabels(ids, tags); + pgen_stats->addColumnValues(0,pgen,mask); + qgen_stats->addColumnValues(0,qgen,mask); + } + + // Find flow parameters for all branch lines + v_vals.clear(); + ids.clear(); + tags.clear(); + mask.clear(); + std::vector id1; + std::vector id2; + std::vector pmin, pmax; + v_vals = pf_app.writeBranchString("flow_str"); + nsize = v_vals.size(); + // Parse branch line endpoints as well as line IDs and values of P and Q for + // base case + for (i=0; i tokens = util.blankTokenizer(v_vals[i]); + if (tokens.size()%8 != 0) { + printf("Incorrect branch power flow listing\n"); + continue; + } + int nline = tokens.size()/8; + for (j=0; jaddRowLabels(id1, id2, tags); + qflow_stats->addRowLabels(id1, id2, tags); + perf_stats->addRowLabels(id1, id2, tags); + pflow_stats->addColumnValues(0,pflow,mask); + qflow_stats->addColumnValues(0,qflow,mask); + perf_stats->addColumnValues(0,perf,mask); + pflow_stats->addRowMinValue(pmin); + qflow_stats->addRowMinValue(pmin); + pflow_stats->addRowMaxValue(pmax); + qflow_stats->addRowMaxValue(pmax); + } + timer->stop(t_store); } - nsize = pflow.size(); - world.max(&nsize,1); -#endif - // Create StatBlock objects for flow parameters and add labels and base case - // values -#ifdef USE_STATBLOCK - gridpack::analysis::StatBlock pflow_stats(world,nsize,ntasks+1); - gridpack::analysis::StatBlock qflow_stats(world,nsize,ntasks+1); - gridpack::analysis::StatBlock perf_stats(world,nsize,ntasks+1); - if (world.rank() == 0) { - pflow_stats.addRowLabels(id1, id2, tags); - qflow_stats.addRowLabels(id1, id2, tags); - perf_stats.addRowLabels(id1, id2, tags); - pflow_stats.addColumnValues(0,pflow,mask); - qflow_stats.addColumnValues(0,qflow,mask); - perf_stats.addColumnValues(0,perf,mask); - pflow_stats.addRowMinValue(pmin); - qflow_stats.addRowMinValue(pmin); - pflow_stats.addRowMaxValue(pmax); - qflow_stats.addRowMaxValue(pmax); - } - timer->stop(t_store); -#endif if (check_Qlim) pf_app.clearQlimViolations(); // Clear any Q limit warnings from base case before starting contingencies gridpack::powerflow::PFBus::clearQlimWarnings(); @@ -744,13 +1835,97 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) // Local contingency results storage for JSON/CSV export std::vector localContingencies; + // Base-case per-element MVA lookup, keyed on (from,to,ckt). + // Populated on rank 0 of each task_comm (only place collectResults ran). + // Used to fill BranchViolation.baseMva/deltaMva during contingency reporting. + std::map baseMvaByKey; + if (outputFormat == "json" || outputFormat == "csv" || + outputFormat == "text") { + for (size_t bi = 0; bi < baseCaseResults.branches.size(); bi++) { + const gridpack::utility::BranchResult &br = baseCaseResults.branches[bi]; + BranchKey k; k.from = br.fromBus; k.to = br.toBus; k.ckt = br.circuitId; + while (!k.ckt.empty() && k.ckt[k.ckt.size()-1] == ' ') k.ckt.resize(k.ckt.size()-1); + double mva = (br.mvaFrom > br.mvaTo) ? br.mvaFrom : br.mvaTo; + baseMvaByKey[k] = mva; + } + } + + // Fill BranchViolation/VoltageViolation arrays from a solved ct result, + // and stream the same rows to _violations..part. + // event_idx of 0 is reserved for base case; contingencies get task_id+1. + auto populateViolations = [&](gridpack::utility::ContingencyResult &ct, + int event_idx) { + const double threshPct = violationSeverityThreshold * 100.0; + // Gate PI on task_comm rank 0 so groupSize>1 doesn't double-count. + bool accumPi = (task_comm.rank() == 0); + for (size_t bi = 0; bi < ct.solution.branches.size(); bi++) { + const gridpack::utility::BranchResult &br = ct.solution.branches[bi]; + if (accumPi) { + double mva_br = (br.mvaFrom > br.mvaTo) ? br.mvaFrom : br.mvaTo; + accumBranchPi(ct.name, mva_br, br.rateSelected); + } + if (br.loadingPercent <= threshPct) continue; + gridpack::utility::BranchViolation v; + v.fromBus = br.fromBus; + v.toBus = br.toBus; + v.circuitId = br.circuitId; + v.mva = (br.mvaFrom > br.mvaTo) ? br.mvaFrom : br.mvaTo; + v.rate = br.rateSelected; + v.loadingPercent = br.loadingPercent; + BranchKey k; k.from = br.fromBus; k.to = br.toBus; k.ckt = br.circuitId; + while (!k.ckt.empty() && k.ckt[k.ckt.size()-1] == ' ') k.ckt.resize(k.ckt.size()-1); + std::map::const_iterator it = baseMvaByKey.find(k); + v.baseMva = (it != baseMvaByKey.end()) ? it->second : 0.0; + v.deltaMva = v.mva - v.baseMva; + v.severity = (br.loadingPercent >= 105.0) ? "critical" : "warning"; + ct.branchViolations.push_back(v); + emitBranchViolation(event_idx, ct.name, v.fromBus, v.toBus, v.circuitId, + v.mva, v.rate, v.loadingPercent, v.baseMva); + } + for (size_t bi = 0; bi < ct.solution.buses.size(); bi++) { + const gridpack::utility::BusResult &b = ct.solution.buses[bi]; + double v_pu = b.voltage; + if (v_pu <= 0.0) continue; // Skip isolated / not-solved buses. + // Voltage PI accrues on every energized bus (textbook form), + // gated on task_comm rank 0 to avoid double-count under groupSize>1. + if (accumPi) accumVoltagePi(ct.name, v_pu); + bool lo = v_pu < Vmin, hi = v_pu > Vmax; + if (!lo && !hi) continue; + gridpack::utility::VoltageViolation vv; + vv.busId = b.busId; + vv.vPu = v_pu; + vv.limitLow = Vmin; + vv.limitHigh = Vmax; + vv.deviationPu = lo ? (v_pu - Vmin) : (v_pu - Vmax); + double dev = std::abs(vv.deviationPu); + vv.severity = (dev >= 0.05) ? "critical" : "warning"; + ct.voltageViolations.push_back(vv); + emitVoltageViolation(event_idx, ct.name, vv.busId, vv.vPu, + vv.limitLow, vv.limitHigh); + } + }; + + // Convergence row recorder; indexes events[task_id]. + auto recordConv = [&](int task_id, const char *status, + const std::string &) { + if (!emitConv) return; + if (task_comm.rank() != 0) return; + ConvRow r; + r.event_idx = task_id + 1; + r.name = events[task_id].p_name; + r.type = (events[task_id].p_type == Branch) ? "branch" : "generator"; + r.cs = pf_app.getConvergence(); + r.status = status; + localConvRows.push_back(r); + }; + // Evaluate contingencies using the task manager int task_id; - char sbuf[128]; + char sbuf[512]; // nextTask returns the same task_id on all processors in task_comm. When the // calculation runs out of task, nextTask will return false. while (taskmgr.nextTask(task_comm, &task_id)) { - printf("Executing task %d on process %d\n",task_id,world.rank()); + if (print_calcs) printf("Executing task %d on process %d\n",task_id,world.rank()); // Trim trailing spaces from contingency name for filename std::string fname = events[task_id].p_name; size_t end = fname.find_last_not_of(' '); @@ -770,7 +1945,7 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) sprintf(sbuf," Line: (from) %d (to) %d (line) \'%s\'\n", events[task_id].p_from[j],events[task_id].p_to[j], events[task_id].p_ckt[j].c_str()); - printf("p[%d] Line: (from) %d (to) %d (line) \'%s\'\n", + if (print_calcs) printf("p[%d] Line: (from) %d (to) %d (line) \'%s\'\n", pf_network->communicator().rank(), events[task_id].p_from[j],events[task_id].p_to[j], events[task_id].p_ckt[j].c_str()); @@ -781,7 +1956,7 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) for (j=0; jcommunicator().rank(), events[task_id].p_busid[j],events[task_id].p_genid[j].c_str()); } @@ -803,10 +1978,6 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) int islandCount = pf_app.getIslandCount(); bool hasLoneBus = pf_app.hasLoneBus(); bool islandDetected = (islandCount > 1); - // Solve power flow equations for this system -#ifdef USE_SUCCESS - contingency_idx.push_back(task_id); -#endif // Skip power flow if contingency setup failed (no valid slack) or islanding detected bool slackCapacityOk = true; // Will be checked after solve bool solveOk = false; @@ -838,52 +2009,52 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) if (!slackCapacityOk) { // Slack generator exceeds Pmax - insufficient generation capacity // This is treated as a failure, similar to divergence -#ifdef USE_SUCCESS - contingency_success.push_back(false); - contingency_violation.push_back(0); - contingency_isolated.push_back(false); -#endif - if (outputFormat != "text") { + if (outputFormat == "json" || outputFormat == "csv") { gridpack::utility::ContingencyResult ctResult; ctResult.name = events[task_id].p_name; ctResult.type = (events[task_id].p_type == Branch) ? "branch" : "generator"; ctResult.hasVoltageViolation = false; ctResult.hasBranchViolation = false; + ctResult.solution.convergence = pf_app.getConvergence(); ctResult.solution.convergence.converged = false; localContingencies.push_back(ctResult); } + recordConv(task_id, "SLACK_OVERLOAD", std::string()); sprintf(sbuf,"\nInsufficient generation capacity for contingency %s\n", events[task_id].p_name.c_str()); if (print_calcs) pf_app.print(sbuf); } else { // Power flow solved and slack within capacity -#ifdef USE_SUCCESS - contingency_success.push_back(true); - contingency_isolated.push_back(hasLoneBus); -#endif // If power flow solution is successful, write out voltages and currents if (print_calcs) pf_app.write(); // Check for violations bool ok1 = pf_app.checkVoltageViolations(); bool ok2 = pf_app.checkLineOverloadViolations(); bool ok = ok1 && ok2; - // Collect results for JSON/CSV export - if (outputFormat != "text") { + // text mode runs the summary path but discards the per-ct struct. + if (outputFormat == "json" || outputFormat == "csv" || + outputFormat == "text") { gridpack::utility::ContingencyResult ctResult; ctResult.name = events[task_id].p_name; ctResult.type = (events[task_id].p_type == Branch) ? "branch" : "generator"; ctResult.hasVoltageViolation = !ok1; ctResult.hasBranchViolation = !ok2; ctResult.solution = pf_app.collectResults(); - localContingencies.push_back(ctResult); + filterResults(ctResult.solution); + populateViolations(ctResult, static_cast(task_id) + 1); + if (outputFormat != "text") localContingencies.push_back(ctResult); + } + if (outputFormat == "csv_flat") { + captureFlatRows(task_id + 1, events[task_id].p_name, true, false); + } + if (outputFormat == "csv_delta") { + captureDeltaRows(task_id + 1, events[task_id], true); } + recordConv(task_id, "OK", std::string()); // Include results of violation checks in output if (ok) { sprintf(sbuf,"\nNo violation for contingency %s\n", events[task_id].p_name.c_str()); -#ifdef USE_SUCCESS - contingency_violation.push_back(1); -#endif } // Report bus voltage violations if (!ok1) { @@ -897,137 +2068,133 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) if (print_calcs) pf_app.writeCABus(); // Report branch overload violations if (!ok2) { - sprintf(sbuf,"\nBranch Violation for contingency %s\n", + // Keep in step with the row format in PFBranch::serialWrite("flow"). + sprintf(sbuf,"\nBranch Violation for contingency %s\n" + " From Bus To Bus CKT P_from Q_from" + " MVA_from P_to Q_to MVA_to" + " Rate Loading%%\n", events[task_id].p_name.c_str()); } else if (!ok) { sprintf(sbuf,"\nNo Branch Violation for contingency %s\n", events[task_id].p_name.c_str()); } -#ifdef USE_SUCCESS - if (!ok1 && !ok2) { - contingency_violation.push_back(4); - } else if (!ok1) { - contingency_violation.push_back(2); - } else if (!ok2) { - contingency_violation.push_back(3); - } -#endif - if (print_calcs) pf_app.print(sbuf); if (print_calcs) pf_app.writeCABranch(); // Get strings of data from power flow calculation and parse them to // extract numerical values. Store these values in vectors and then // add them to StatBlock objects -#ifdef USE_STATBLOCK - timer->start(t_store); - vmag.clear(); - vang.clear(); - mask.clear(); - mag_mask.clear(); - v_vals.clear(); - v_vals = pf_app.writeBusString("vr_str"); - nsize = v_vals.size(); - for (i=0; i tokens = util.blankTokenizer(v_vals[i]); - int not_isolated = atoi(tokens[3].c_str()); - if (not_isolated == 1) { - vmag.push_back(atof(tokens[2].c_str())); - if (atoi(tokens[4].c_str()) != 0) { - mag_mask.push_back(2); - } else { - mag_mask.push_back(1); + if (write_stats) { + timer->start(t_store); + vmag.clear(); + vang.clear(); + mask.clear(); + mag_mask.clear(); + v_vals.clear(); + v_vals = pf_app.writeBusString("vr_str"); + nsize = v_vals.size(); + for (i=0; i tokens = util.blankTokenizer(v_vals[i]); + if (!busMonitored(atoi(tokens[0].c_str()))) continue; + int not_isolated = atoi(tokens[3].c_str()); + if (not_isolated == 1) { + vmag.push_back(atof(tokens[2].c_str())); + if (atoi(tokens[4].c_str()) != 0) { + mag_mask.push_back(2); + } else { + mag_mask.push_back(1); + } } - } - vang.push_back(atof(tokens[1].c_str())); - mask.push_back(1); - } -#endif -#ifdef USE_STATBLOCK - if (task_comm.rank() == 0) { - vmag_stats.addColumnValues(task_id+1,vmag,mag_mask); - vang_stats.addColumnValues(task_id+1,vang,mask); - } -#endif -#ifdef USE_STATBLOCK - pgen.clear(); - qgen.clear(); - mask.clear(); - v_vals.clear(); - v_vals = pf_app.writeBusString("power"); - nsize = v_vals.size(); - for (i=0; i tokens = util.blankTokenizer(v_vals[i]); - if (tokens.size()%4 != 0) { - printf("Incorrect generator listing\n"); - continue; - } - int ngen = tokens.size()/4; - for (j=0; j tokens = util.blankTokenizer(v_vals[i]); - if (tokens.size()%8 != 0) { - printf("Incorrect branch power flow listing\n"); - continue; + if (task_comm.rank() == 0) { + vmag_stats->addColumnValues(task_id+1,vmag,mag_mask); + vang_stats->addColumnValues(task_id+1,vang,mask); } - int nline = tokens.size()/8; - for (j=0; j tokens = util.blankTokenizer(v_vals[i]); + if (tokens.size()%4 != 0) { + printf("Incorrect generator listing\n"); + continue; + } + int ngen = tokens.size()/4; + for (j=0; jaddColumnValues(task_id+1,pgen,mask); + qgen_stats->addColumnValues(task_id+1,qgen,mask); + } + pflow.clear(); + qflow.clear(); + perf.clear(); + mask.clear(); + v_vals.clear(); + v_vals = pf_app.writeBranchString("flow_str"); + nsize = v_vals.size(); + for (i=0; i tokens = util.blankTokenizer(v_vals[i]); + if (tokens.size()%8 != 0) { + printf("Incorrect branch power flow listing\n"); + continue; + } + int nline = tokens.size()/8; + for (j=0; jaddColumnValues(task_id+1,pflow,mask); + qflow_stats->addColumnValues(task_id+1,qflow,mask); + perf_stats->addColumnValues(task_id+1,perf,mask); + } + timer->stop(t_store); } -#endif -#ifdef USE_STATBLOCK - if (task_comm.rank() == 0) { - pflow_stats.addColumnValues(task_id+1,pflow,mask); - qflow_stats.addColumnValues(task_id+1,qflow,mask); - perf_stats.addColumnValues(task_id+1,perf,mask); - } - timer->stop(t_store); -#endif // Note: clearQlimViolations() moved after unSetContingency() below } // end slackCapacityOk block } else { -#ifdef USE_SUCCESS - contingency_success.push_back(false); - contingency_violation.push_back(0); - contingency_isolated.push_back(false); -#endif - if (outputFormat != "text") { + if (outputFormat == "json" || outputFormat == "csv") { gridpack::utility::ContingencyResult ctResult; ctResult.name = events[task_id].p_name; ctResult.type = (events[task_id].p_type == Branch) ? "branch" : "generator"; ctResult.hasVoltageViolation = false; ctResult.hasBranchViolation = false; + ctResult.solution.convergence = pf_app.getConvergence(); ctResult.solution.convergence.converged = false; localContingencies.push_back(ctResult); } + { + const char *st; + if (islandDetected) { + st = "ISLANDED"; + } else if (!contingencyFound) { + st = "NO_SLACK"; + } else { + st = "DIVERGED"; + } + recordConv(task_id, st, std::string()); + } if (islandDetected) { sprintf(sbuf,"\nIslanding detected for contingency %s (%d islands)\n", events[task_id].p_name.c_str(), islandCount); @@ -1041,90 +2208,85 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) if (print_calcs) pf_app.print(sbuf); // Add dummy values to StatBlock object. Mask value is set to 0 for all // network elements to indicate calculation failure -#ifdef USE_STATBLOCK - timer->start(t_store); - vmag.clear(); - vang.clear(); - mask.clear(); - mag_mask.clear(); - v_vals.clear(); - v_vals = pf_app.writeBusString("vfail_str"); - nsize = v_vals.size(); - for (i=0; i tokens = util.blankTokenizer(v_vals[i]); - int not_isolated = atoi(tokens[3].c_str()); - if (not_isolated == 1) { - vmag.push_back(0.0); - mag_mask.push_back(0); + if (write_stats) { + timer->start(t_store); + vmag.clear(); + vang.clear(); + mask.clear(); + mag_mask.clear(); + v_vals.clear(); + v_vals = pf_app.writeBusString("vfail_str"); + nsize = v_vals.size(); + for (i=0; i tokens = util.blankTokenizer(v_vals[i]); + if (!busMonitored(atoi(tokens[0].c_str()))) continue; + int not_isolated = atoi(tokens[3].c_str()); + if (not_isolated == 1) { + vmag.push_back(0.0); + mag_mask.push_back(0); + } + vang.push_back(0.0); + mask.push_back(0); } - vang.push_back(0.0); - mask.push_back(0); - } -#endif -#ifdef USE_STATBLOCK - if (task_comm.rank() == 0) { - vmag_stats.addColumnValues(task_id+1,vmag,mag_mask); - vang_stats.addColumnValues(task_id+1,vang,mask); - } -#endif -#ifdef USE_STATBLOCK - pgen.clear(); - qgen.clear(); - mask.clear(); - v_vals.clear(); - v_vals = pf_app.writeBusString("pfail_str"); - nsize = v_vals.size(); - for (i=0; i tokens = util.blankTokenizer(v_vals[i]); - if (tokens.size()%4 != 0) { - printf("Incorrect generator listing\n"); - continue; + if (task_comm.rank() == 0) { + vmag_stats->addColumnValues(task_id+1,vmag,mag_mask); + vang_stats->addColumnValues(task_id+1,vang,mask); } - int ngen = tokens.size()/4; - for (j=0; j tokens = util.blankTokenizer(v_vals[i]); + if (tokens.size()%4 != 0) { + printf("Incorrect generator listing\n"); + continue; + } + int ngen = tokens.size()/4; + for (j=0; j tokens = util.blankTokenizer(v_vals[i]); - if (tokens.size()%8 != 0) { - printf("Incorrect branch power flow listing\n"); - continue; + if (task_comm.rank() == 0) { + pgen_stats->addColumnValues(task_id+1,pgen,mask); + qgen_stats->addColumnValues(task_id+1,qgen,mask); } - int nline = tokens.size()/8; - for (j=0; j tokens = util.blankTokenizer(v_vals[i]); + if (tokens.size()%8 != 0) { + printf("Incorrect branch power flow listing\n"); + continue; + } + int nline = tokens.size()/8; + for (j=0; jaddColumnValues(task_id+1,pflow,mask); + qflow_stats->addColumnValues(task_id+1,qflow,mask); + perf_stats->addColumnValues(task_id+1,perf,mask); + } + timer->stop(t_store); } -#endif -#ifdef USE_STATBLOCK - if (task_comm.rank() == 0) { - pflow_stats.addColumnValues(task_id+1,pflow,mask); - qflow_stats.addColumnValues(task_id+1,qflow,mask); - perf_stats.addColumnValues(task_id+1,perf,mask); - } - timer->stop(t_store); -#endif } // Return network to its original base case state pf_app.unSetContingency(events[task_id]); @@ -1137,58 +2299,506 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) // Close output file for this contingency if (print_calcs) pf_app.close(); } - // Print statistics from task manager describing the number of tasks performed - // per processor - taskmgr.printStats(); + // csv_flat / csv_delta: each rank streamed rows to its .part file during + // the loop. Close, sync, then world rank 0 writes header + concatenates. + if (outputFormat == "csv_flat") { + if (flatPart.is_open()) flatPart.close(); + } + if (outputFormat == "csv_delta") { + if (deltaPart.is_open()) deltaPart.close(); + } + if (wantBusSidecar) { + world.sync(); + if (world.rank() == 0) { + const size_t BUFSZ = 1 << 20; + std::vector buf(BUFSZ); + + auto concatParts = [&](const char *suffix, const char *header, + const char *tag, const char *outName) { + std::string outFile = outputFile + outName; + std::ofstream fout(outFile.c_str(), + std::ios::out | std::ios::trunc | std::ios::binary); + fout << header; + size_t rows = 0; + for (int p = 0; p < world.size(); p++) { + std::ostringstream oss; + oss << outputFile << suffix << p << ".part"; + std::string part = oss.str(); + std::ifstream fin(part.c_str(), std::ios::in | std::ios::binary); + if (!fin) continue; + while (fin) { + fin.read(&buf[0], BUFSZ); + std::streamsize got = fin.gcount(); + if (got > 0) { + fout.write(&buf[0], got); + for (std::streamsize k = 0; k < got; k++) { + if (buf[k] == '\n') rows++; + } + } + } + fin.close(); + std::remove(part.c_str()); + } + fout.close(); + printf("[%s] wrote %zu rows to %s\n", tag, rows, outFile.c_str()); + }; + + if (outputFormat == "csv_flat") { + concatParts("_flat.", + "event_idx,contingency,from_bus,to_bus,circuit_id," + "p_from_mw,q_from_mvar,mva_from,rate_mva,loading_percent," + "viol,v_from_pu,v_to_pu,ang_from_deg,ang_to_deg\n", + "csv_flat", + "_flat.csv"); + } + if (outputFormat == "csv_delta") { + concatParts("_delta.", + "event_idx,contingency,type,from_bus,to_bus,ckt," + "base_kv_from,base_kv_to,area_from,area_to,base_rate_mva,cont_rate_mva," + "base_p_mw,cont_p_mw,base_q_mvar,cont_q_mvar," + "base_mva,cont_mva,base_loading_pct,cont_loading_pct," + "v_from_base,v_from_cont,v_to_base,v_to_cont," + "ang_from_base,ang_from_cont,ang_to_base,ang_to_cont," + // ",cont_event_facility" here if re-enabling the column + "d_v_base,d_v_cont,d_angle_base,d_angle_cont\n", + "csv_delta", + "_delta.csv"); + } - // Gather stats on successful contingency calculations -#ifdef USE_SUCCESS - if (task_comm.rank() == 0) { - ca_success.addElements(contingency_idx, contingency_success); - ca_violation.addElements(contingency_idx, contingency_violation); - ca_isolated.addElements(contingency_idx, contingency_isolated); + // Bus metadata sidecar (deduped by bus_id, first writer wins). + std::string busFile = outputFile + "_buses.csv"; + std::ofstream bout(busFile.c_str(), + std::ios::out | std::ios::trunc | std::ios::binary); + bout << "bus_id,bus_name,base_kv,area,zone,owner," + "area_name,zone_name,owner_name\n"; + std::set seen_bus; + size_t bus_rows = 0; + for (int p = 0; p < world.size(); p++) { + std::ostringstream oss; + oss << outputFile << "_buses." << p << ".part"; + std::string part = oss.str(); + std::ifstream fin(part.c_str()); + if (!fin) continue; + std::string line; + while (std::getline(fin, line)) { + if (line.empty()) continue; + size_t comma = line.find(','); + if (comma == std::string::npos) continue; + int bus_id = std::atoi(line.substr(0, comma).c_str()); + if (seen_bus.insert(bus_id).second) { + bout << line << "\n"; + bus_rows++; + } + } + fin.close(); + std::remove(part.c_str()); + } + bout.close(); + printf("[buses] wrote %zu rows to %s\n", bus_rows, busFile.c_str()); + } } - ca_success.upload(); - ca_violation.upload(); - ca_isolated.upload(); - // All processes call getData to ensure GA progress (NGA_Gather requires - // remote process participation for one-sided communication). - contingency_idx.clear(); - contingency_success.clear(); - contingency_violation.clear(); - contingency_isolated.clear(); - for (i=0; i.part into _violations.csv. + // Runs for every outputFormat; ranks that emitted zero rows simply have no + // .part file to include. + if (violPart.is_open()) violPart.close(); + world.sync(); if (world.rank() == 0) { - std::ofstream fout; - fout.open("success.txt"); - for (i=0; i buf(BUFSZ); + std::string outFile = outputFile + "_violations.csv"; + std::ofstream fout(outFile.c_str(), + std::ios::out | std::ios::trunc | std::ios::binary); + fout << "event_idx,contingency,type,element,mva_or_vpu,rate_or_limit," + "loading_percent,base_mva,delta,severity\n"; + size_t rows = 0; + for (int p = 0; p < world.size(); p++) { + std::ostringstream oss; + oss << outputFile << "_violations." << p << ".part"; + std::string part = oss.str(); + std::ifstream fin(part.c_str(), std::ios::in | std::ios::binary); + if (!fin) continue; + while (fin) { + fin.read(&buf[0], BUFSZ); + std::streamsize got = fin.gcount(); + if (got > 0) { + fout.write(&buf[0], got); + for (std::streamsize k = 0; k < got; k++) { + if (buf[k] == '\n') rows++; + } } - if (contingency_isolated[i]) { - fout << " warning: isolated"; + } + fin.close(); + std::remove(part.c_str()); + } + fout.close(); + printf("[violations] wrote %zu rows to %s\n", rows, outFile.c_str()); + } + + // Aggregate skip count across ranks for diagnostics. + if (outputFormat == "csv_delta") { + long localSkip = static_cast(deltaSkipCount); + long totalSkip = localSkip; + world.sum(&totalSkip, 1); + if (world.rank() == 0 && totalSkip > 0) { + printf("[csv_delta] %ld branch rows had no base-cache match\n", + totalSkip); + } + } + + // Aggregate per-run summary across all ranks. Emit _summary.json + // on rank 0. Counters and worst-of values follow the same shape commercial + // tools use in their contingency reports. + { + // 0: total_ct 1: converged 2: cts_with_branch_viol 3: cts_with_voltage_viol + // 4: islanded 5: no_slack 6: diverged 7: slack_overload + long localCounters[8] = { 0 }; + if (!localContingencies.empty()) { + // json/csv paths: authoritative per-ct list. + for (size_t ci = 0; ci < localContingencies.size(); ci++) { + const gridpack::utility::ContingencyResult &ct = localContingencies[ci]; + localCounters[0] += 1; + if (ct.solution.convergence.converged) localCounters[1] += 1; + } + } else { + // text/csv_flat/csv_delta: count from convergence rows. + localCounters[0] = static_cast(localConvRows.size()); + long conv = 0; + for (size_t i = 0; i < localConvRows.size(); i++) { + if (localConvRows[i].status == "OK") conv++; + } + localCounters[1] = conv; + } + localCounters[2] = static_cast(ctsWithBranchViol.size()); + localCounters[3] = static_cast(ctsWithVoltageViol.size()); + // Per-status breakdown from convergence rows (populated in every mode). + for (size_t i = 0; i < localConvRows.size(); i++) { + const std::string &st = localConvRows[i].status; + if (st == "ISLANDED") localCounters[4] += 1; + else if (st == "NO_SLACK") localCounters[5] += 1; + else if (st == "DIVERGED") localCounters[6] += 1; + else if (st == "SLACK_OVERLOAD") localCounters[7] += 1; + } + long totalCounters[8] = { 0 }; + for (int i = 0; i < 8; i++) totalCounters[i] = localCounters[i]; + world.sum(&totalCounters[0], 8); + long localViolRows = static_cast(violRowCount); + long totalViolRows = localViolRows; + world.sum(&totalViolRows, 1); + + struct WorstBranchWire { + double loading_pct; + int from, to; + char ckt[4]; + char ct_name[32]; + }; + struct WorstVoltageWire { + double v_pu; + double dev_pu; + int bus_id; + char ct_name[32]; + }; + WorstBranchWire wbLo; + wbLo.loading_pct = worstBr.loading_pct; + wbLo.from = worstBr.from; wbLo.to = worstBr.to; + std::strncpy(wbLo.ckt, worstBr.ckt.c_str(), 3); wbLo.ckt[3] = '\0'; + std::strncpy(wbLo.ct_name, worstBr.ct_name.c_str(), 31); wbLo.ct_name[31] = '\0'; + WorstVoltageWire wvLo; + wvLo.v_pu = worstVLo.v_pu; wvLo.dev_pu = worstVLo.dev_pu; + wvLo.bus_id = worstVLo.bus_id; + std::strncpy(wvLo.ct_name, worstVLo.ct_name.c_str(), 31); wvLo.ct_name[31] = '\0'; + WorstVoltageWire wvHi; + wvHi.v_pu = worstVHi.v_pu; wvHi.dev_pu = worstVHi.dev_pu; + wvHi.bus_id = worstVHi.bus_id; + std::strncpy(wvHi.ct_name, worstVHi.ct_name.c_str(), 31); wvHi.ct_name[31] = '\0'; + MPI_Comm mpi_comm = static_cast(world); + if (world.rank() == 0) { + for (int p = 1; p < world.size(); p++) { + WorstBranchWire otherB; + WorstVoltageWire otherLo, otherHi; + MPI_Recv(&otherB, sizeof(otherB), MPI_BYTE, p, 30, mpi_comm, MPI_STATUS_IGNORE); + MPI_Recv(&otherLo, sizeof(otherLo), MPI_BYTE, p, 31, mpi_comm, MPI_STATUS_IGNORE); + MPI_Recv(&otherHi, sizeof(otherHi), MPI_BYTE, p, 32, mpi_comm, MPI_STATUS_IGNORE); + if (otherB.loading_pct > wbLo.loading_pct) wbLo = otherB; + if (otherLo.dev_pu < wvLo.dev_pu) wvLo = otherLo; + if (otherHi.dev_pu > wvHi.dev_pu) wvHi = otherHi; + } + } else { + MPI_Send(&wbLo, sizeof(wbLo), MPI_BYTE, 0, 30, mpi_comm); + MPI_Send(&wvLo, sizeof(wvLo), MPI_BYTE, 0, 31, mpi_comm); + MPI_Send(&wvHi, sizeof(wvHi), MPI_BYTE, 0, 32, mpi_comm); + } + + // Gather per-ct PI/VPI/Vdev/rosters + per-ct worst-single-element data. + // Blob line: name pi vpi vdev worstLoad worstVdevAbs worstVpu br_flag v_flag + std::map aggPi; + std::map aggVpi; + std::map aggVdev; + std::map aggWorstLoading; + std::map aggWorstVdev; // unsigned max |v-limit| + std::map aggWorstVpu; // v_pu that produced aggWorstVdev + std::set aggBranchViol; + std::set aggVoltageViol; + { + std::ostringstream localOut; + localOut << std::setprecision(10); + std::set names; + for (std::map::const_iterator it = ctPi.begin(); + it != ctPi.end(); ++it) names.insert(it->first); + for (std::map::const_iterator it = ctVpi.begin(); + it != ctVpi.end(); ++it) names.insert(it->first); + for (std::map::const_iterator it = ctVdev.begin(); + it != ctVdev.end(); ++it) names.insert(it->first); + for (std::set::const_iterator it = ctsWithBranchViol.begin(); + it != ctsWithBranchViol.end(); ++it) names.insert(*it); + for (std::set::const_iterator it = ctsWithVoltageViol.begin(); + it != ctsWithVoltageViol.end(); ++it) names.insert(*it); + auto lookup = [](const std::map &m, + const std::string &k) -> double { + std::map::const_iterator it = m.find(k); + return (it == m.end()) ? 0.0 : it->second; + }; + for (std::set::const_iterator it = names.begin(); + it != names.end(); ++it) { + const std::string &nm = *it; + double pi = lookup(ctPi, nm); + double vpi = lookup(ctVpi, nm); + double vdev = lookup(ctVdev, nm); + double wL = lookup(ctWorstLoading, nm); + double wDabs = lookup(ctWorstVdev, nm); + double wV = lookup(ctWorstVpu, nm); + int brFlag = ctsWithBranchViol.count(nm) ? 1 : 0; + int vFlag = ctsWithVoltageViol.count(nm) ? 1 : 0; + localOut << nm << '\t' << pi << '\t' << vpi << '\t' << vdev << '\t' + << wL << '\t' << wDabs << '\t' << wV << '\t' + << brFlag << '\t' << vFlag << '\n'; + } + std::string localBlob = localOut.str(); + auto splitTabs = [](const std::string &line, + std::vector &out) { + out.clear(); + size_t pos = 0; + while (pos <= line.size()) { + size_t t = line.find('\t', pos); + if (t == std::string::npos) { + out.push_back(line.substr(pos)); + break; + } + out.push_back(line.substr(pos, t - pos)); + pos = t + 1; + } + }; + auto absorb = [&](const std::string &blob) { + size_t pos = 0; + std::vector fields; + while (pos < blob.size()) { + size_t eol = blob.find('\n', pos); + if (eol == std::string::npos) break; + std::string line = blob.substr(pos, eol - pos); + pos = eol + 1; + splitTabs(line, fields); + if (fields.size() < 9) continue; + const std::string &nm = fields[0]; + double pi = std::atof(fields[1].c_str()); + double vpi = std::atof(fields[2].c_str()); + double vdev = std::atof(fields[3].c_str()); + double wL = std::atof(fields[4].c_str()); + double wDabs = std::atof(fields[5].c_str()); + double wV = std::atof(fields[6].c_str()); + int brF = std::atoi(fields[7].c_str()); + int vF = std::atoi(fields[8].c_str()); + if (pi != 0.0) aggPi[nm] += pi; + if (vpi != 0.0) aggVpi[nm] += vpi; + if (vdev != 0.0) aggVdev[nm] += vdev; + if (wL > aggWorstLoading[nm]) aggWorstLoading[nm] = wL; + if (wDabs > aggWorstVdev[nm]) { aggWorstVdev[nm] = wDabs; aggWorstVpu[nm] = wV; } + if (brF) aggBranchViol.insert(nm); + if (vF) aggVoltageViol.insert(nm); + } + }; + if (world.rank() == 0) { + absorb(localBlob); + for (int p = 1; p < world.size(); p++) { + int len = 0; + MPI_Recv(&len, 1, MPI_INT, p, 33, mpi_comm, MPI_STATUS_IGNORE); + std::string remote; + remote.resize(len); + if (len > 0) { + MPI_Recv(&remote[0], len, MPI_CHAR, p, 34, mpi_comm, + MPI_STATUS_IGNORE); + } + absorb(remote); } - fout << std::endl; } else { - fout << "contingency: " << i+1 << " success: false" << std::endl; + int len = static_cast(localBlob.size()); + MPI_Send(&len, 1, MPI_INT, 0, 33, mpi_comm); + if (len > 0) { + MPI_Send(const_cast(localBlob.c_str()), len, MPI_CHAR, 0, 34, + mpi_comm); + } } } - fout.close(); + + if (world.rank() == 0) { + std::string sumFile = outputFile + "_summary.json"; + std::ofstream sout(sumFile.c_str()); + sout << std::fixed; + sout << "{\n"; + sout << " \"total_contingencies\": " << totalCounters[0] << ",\n"; + sout << " \"converged\": " << totalCounters[1] << ",\n"; + sout << " \"diverged\": " << (totalCounters[0] - totalCounters[1]) << ",\n"; + // Per-status breakdown of the diverged bucket. Sum equals `diverged`. + sout << " \"islanded\": " << totalCounters[4] << ",\n"; + sout << " \"no_slack\": " << totalCounters[5] << ",\n"; + sout << " \"solver_diverged\": " << totalCounters[6] << ",\n"; + sout << " \"slack_overload\": " << totalCounters[7] << ",\n"; + sout << " \"with_branch_violation\": " << totalCounters[2] << ",\n"; + sout << " \"with_voltage_violation\": " << totalCounters[3] << ",\n"; + sout << " \"worst_loading\": "; + if (wbLo.loading_pct > 0.0) { + sout << "{\"contingency\": \"" << wbLo.ct_name << "\"" + << ", \"from_bus\": " << wbLo.from + << ", \"to_bus\": " << wbLo.to + << ", \"circuit_id\": \"" << wbLo.ckt << "\"" + << ", \"loading_percent\": " << std::setprecision(2) << wbLo.loading_pct + << "},\n"; + } else { + sout << "null,\n"; + } + sout << " \"worst_voltage_low\": "; + if (wvLo.dev_pu < 0.0) { + sout << "{\"contingency\": \"" << wvLo.ct_name << "\"" + << ", \"bus_id\": " << wvLo.bus_id + << ", \"v_pu\": " << std::setprecision(6) << wvLo.v_pu + << ", \"deviation_pu\": " << std::setprecision(6) << wvLo.dev_pu + << "},\n"; + } else { + sout << "null,\n"; + } + sout << " \"worst_voltage_high\": "; + if (wvHi.dev_pu > 0.0) { + sout << "{\"contingency\": \"" << wvHi.ct_name << "\"" + << ", \"bus_id\": " << wvHi.bus_id + << ", \"v_pu\": " << std::setprecision(6) << wvHi.v_pu + << ", \"deviation_pu\": " << std::setprecision(6) << wvHi.dev_pu + << "},\n"; + } else { + sout << "null,\n"; + } + auto emitNameArray = [&](const char *field, + const std::set &names) { + sout << " \"" << field << "\": ["; + int emitted = 0; + for (std::set::const_iterator it = names.begin(); + it != names.end(); ++it, ++emitted) { + if (emitted) sout << ", "; + sout << "\"" << *it << "\""; + } + sout << "],\n"; + }; + emitNameArray("contingencies_with_branch_violation", aggBranchViol); + emitNameArray("contingencies_with_voltage_violation", aggVoltageViol); + // top_severe_contingencies: Group A (any violation) first, sorted by + // worst-single-element severity; Group B (no violations) after, + // sorted by composite_pi. Combined list capped at topN. + { + std::set ctSet; + for (std::map::const_iterator it = aggPi.begin(); + it != aggPi.end(); ++it) ctSet.insert(it->first); + for (std::map::const_iterator it = aggVpi.begin(); + it != aggVpi.end(); ++it) ctSet.insert(it->first); + for (std::set::const_iterator it = aggBranchViol.begin(); + it != aggBranchViol.end(); ++it) ctSet.insert(*it); + for (std::set::const_iterator it = aggVoltageViol.begin(); + it != aggVoltageViol.end(); ++it) ctSet.insert(*it); + auto agg_get = [](const std::map &m, + const std::string &k) -> double { + std::map::const_iterator it = m.find(k); + return (it == m.end()) ? 0.0 : it->second; + }; + struct SevRow { + std::string name; + bool violated; + double sortKey; // Group A: worst-single severity; Group B: composite_pi + double composite, branchPi, voltagePi; + double worstLoading; + double worstVpu, worstVdev; // worstVdev is unsigned |v - limit| + }; + std::vector groupA, groupB; + for (std::set::const_iterator it = ctSet.begin(); + it != ctSet.end(); ++it) { + SevRow r; + r.name = *it; + r.branchPi = agg_get(aggPi, *it); + r.voltagePi = agg_get(aggVpi, *it); + r.composite = piBranchWeight * r.branchPi + piVoltageWeight * r.voltagePi; + r.worstLoading = agg_get(aggWorstLoading, *it); + r.worstVdev = agg_get(aggWorstVdev, *it); + r.worstVpu = agg_get(aggWorstVpu, *it); + bool hasBr = aggBranchViol.count(*it) > 0; + bool hasV = aggVoltageViol.count(*it) > 0; + r.violated = hasBr || hasV; + if (r.violated) { + double s = 0.0; + if (r.worstLoading > 100.0) s = std::max(s, r.worstLoading - 100.0); + if (r.worstVdev > 0.0) s = std::max(s, r.worstVdev * 1000.0); + r.sortKey = s; + groupA.push_back(r); + } else { + r.sortKey = r.composite; + if (r.sortKey > 0.0) groupB.push_back(r); + } + } + auto sevCmp = [](const SevRow &a, const SevRow &b) { + if (a.sortKey != b.sortKey) return a.sortKey > b.sortKey; + return a.name < b.name; + }; + std::sort(groupA.begin(), groupA.end(), sevCmp); + std::sort(groupB.begin(), groupB.end(), sevCmp); + auto emitRow = [&](const SevRow &r, bool first) { + if (!first) sout << ",\n "; + else sout << "\n "; + bool hasBr = aggBranchViol.count(r.name) > 0; + bool hasV = r.worstVdev > 0.0; + sout << "{\"contingency\": \"" << r.name << "\"" + << ", \"composite_pi\": " << std::setprecision(6) << r.composite + << ", \"worst_branch_loading_percent\": "; + if (hasBr) sout << std::setprecision(2) << r.worstLoading; + else sout << "null"; + sout << ", \"worst_voltage_pu\": "; + if (hasV) sout << std::setprecision(6) << r.worstVpu; + else sout << "null"; + sout << ", \"worst_voltage_deviation_pu\": "; + if (hasV) sout << std::setprecision(6) << r.worstVdev; + else sout << "null"; + sout << ", \"has_branch_violation\": " + << (hasBr ? "true" : "false") + << ", \"has_voltage_violation\": " + << (aggVoltageViol.count(r.name) ? "true" : "false") + << "}"; + }; + sout << " \"top_severe_contingencies\": ["; + int emitted = 0; + for (size_t i = 0; i < groupA.size() && emitted < topN; ++i, ++emitted) { + emitRow(groupA[i], emitted == 0); + } + for (size_t i = 0; i < groupB.size() && emitted < topN; ++i, ++emitted) { + emitRow(groupB[i], emitted == 0); + } + sout << (emitted ? "\n ]\n" : "]\n"); + } + sout << "}\n"; + sout.close(); + printf("[summary] wrote %s (%ld contingencies, %ld converged, " + "%ld with branch violations, %ld with voltage violations)\n", + sumFile.c_str(), + totalCounters[0], totalCounters[1], + totalCounters[2], totalCounters[3]); + } } -#endif + + // Print statistics from task manager describing the number of tasks performed + // per processor + taskmgr.printStats(); // Sync GA before MPI collectives to flush any pending one-sided operations world.sync(); @@ -1252,12 +2862,13 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) } if (outputFormat == "csv") { - // Each process serializes its contingency CSV data into 4 strings - // (buses, branches, generators, convergence rows without headers) - std::ostringstream localBus, localBranch, localGen, localConv; + // Convergence is emitted by the universal sidecar block below. + std::ostringstream localBus, localBranch, localGen; localBus << std::fixed; localBranch << std::fixed; localGen << std::fixed; + // Column layout must match the headers written by + // ResultsExporter::writePFCSV for the base case. for (size_t ci = 0; ci < localContingencies.size(); ci++) { const gridpack::utility::ContingencyResult& ct = localContingencies[ci]; const gridpack::utility::PowerFlowResults& r = ct.solution; @@ -1291,6 +2902,7 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) << std::setprecision(4) << br.mvaFrom << "," << std::setprecision(4) << br.mvaTo << "," << std::setprecision(4) << br.rateA << "," + << std::setprecision(4) << br.rateSelected << "," << std::setprecision(2) << br.loadingPercent << "\n"; } for (size_t gi = 0; gi < r.generators.size(); gi++) { @@ -1304,33 +2916,22 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) << std::setprecision(6) << g.voltageSetpoint << "," << g.status << "\n"; } - localConv << ct.name << "," - << (r.convergence.converged ? "true" : "false") << "," - << r.convergence.iterations << "," - << std::scientific << r.convergence.finalTolerance << "," - << std::fixed - << r.convergence.finalMismatch.maxPBus << "," - << std::setprecision(4) << r.convergence.finalMismatch.maxPMismatch << "," - << r.convergence.finalMismatch.maxQBus << "," - << std::setprecision(4) << r.convergence.finalMismatch.maxQMismatch << "\n"; } // Gather all CSV fragments on rank 0 using point-to-point send/recv MPI_Comm mpi_comm = static_cast(world); std::vector allBus(world.size()), allBranch(world.size()); - std::vector allGen(world.size()), allConv(world.size()); + std::vector allGen(world.size()); allBus[0] = localBus.str(); allBranch[0] = localBranch.str(); allGen[0] = localGen.str(); - allConv[0] = localConv.str(); if (world.rank() == 0) { for (int p = 1; p < world.size(); p++) { - int lens[4]; - MPI_Recv(lens, 4, MPI_INT, p, 0, mpi_comm, MPI_STATUS_IGNORE); + int lens[3]; + MPI_Recv(lens, 3, MPI_INT, p, 0, mpi_comm, MPI_STATUS_IGNORE); allBus[p].resize(lens[0]); allBranch[p].resize(lens[1]); allGen[p].resize(lens[2]); - allConv[p].resize(lens[3]); if (lens[0] > 0) MPI_Recv(&allBus[p][0], lens[0], MPI_CHAR, p, 1, mpi_comm, MPI_STATUS_IGNORE); @@ -1340,16 +2941,12 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) if (lens[2] > 0) MPI_Recv(&allGen[p][0], lens[2], MPI_CHAR, p, 3, mpi_comm, MPI_STATUS_IGNORE); - if (lens[3] > 0) - MPI_Recv(&allConv[p][0], lens[3], MPI_CHAR, p, 4, mpi_comm, - MPI_STATUS_IGNORE); } } else { std::string sBus = localBus.str(), sBranch = localBranch.str(); - std::string sGen = localGen.str(), sConv = localConv.str(); - int lens[4] = {(int)sBus.size(), (int)sBranch.size(), - (int)sGen.size(), (int)sConv.size()}; - MPI_Send(lens, 4, MPI_INT, 0, 0, mpi_comm); + std::string sGen = localGen.str(); + int lens[3] = {(int)sBus.size(), (int)sBranch.size(), (int)sGen.size()}; + MPI_Send(lens, 3, MPI_INT, 0, 0, mpi_comm); if (lens[0] > 0) MPI_Send(const_cast(sBus.c_str()), lens[0], MPI_CHAR, 0, 1, mpi_comm); @@ -1359,9 +2956,6 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) if (lens[2] > 0) MPI_Send(const_cast(sGen.c_str()), lens[2], MPI_CHAR, 0, 3, mpi_comm); - if (lens[3] > 0) - MPI_Send(const_cast(sConv.c_str()), lens[3], MPI_CHAR, 0, 4, - mpi_comm); } // Rank 0 writes the CSV files @@ -1382,52 +2976,119 @@ void gridpack::contingency_analysis::CADriver::execute(int argc, char** argv) std::ofstream out((outputFile + "_generators.csv").c_str(), std::ios::app); for (size_t p = 0; p < allGen.size(); p++) out << allGen[p]; } - { - std::ofstream out((outputFile + "_convergence.csv").c_str(), std::ios::app); - for (size_t p = 0; p < allConv.size(); p++) out << allConv[p]; - } -#ifdef USE_SUCCESS - // Write summary CSV - std::string summaryFile = outputFile + "_summary.csv"; - std::ofstream sout(summaryFile.c_str()); - sout << "contingency,type,converged,has_voltage_violation,has_branch_violation\n"; - for (int ci = 0; ci < ntasks; ci++) { - bool converged = (contingency_violation[ci] > 0); - sout << events[ci].p_name << "," - << (events[ci].p_type == Branch ? "branch" : "generator") << "," - << (converged ? "true" : "false") << "," - << ((contingency_violation[ci] == 2 || contingency_violation[ci] == 4) - ? "true" : "false") << "," - << ((contingency_violation[ci] == 3 || contingency_violation[ci] == 4) - ? "true" : "false") << "\n"; + } + } + + // Universal convergence sidecar: gather, sort by event_idx, write. + if (emitConv) { + auto formatRow = [](std::ostringstream &os, const ConvRow &r) { + // Derive converged from status so ISLANDED/NO_SLACK cases -- where solve() + // was never entered and pf_app.getConvergence() returns the previous + // case's stale value -- read as false, matching _summary.json's + // diverged=total-converged accounting. + os << r.event_idx << "," + << r.name << "," + << r.type << "," + << ((r.status == "OK") ? "true" : "false") << "," + << r.cs.iterations << "," + << std::scientific << r.cs.finalTolerance << "," + << std::fixed + << r.cs.finalMismatch.maxPBus << "," + << std::setprecision(4) << r.cs.finalMismatch.maxPMismatch << "," + << r.cs.finalMismatch.maxQBus << "," + << std::setprecision(4) << r.cs.finalMismatch.maxQMismatch << "," + << r.status << "\n"; + }; + + std::vector idx; + std::ostringstream localStream; + localStream << std::fixed; + std::vector localOffsets; + localOffsets.reserve(localConvRows.size() + 1); + for (size_t i = 0; i < localConvRows.size(); i++) { + localOffsets.push_back(static_cast(localStream.tellp())); + formatRow(localStream, localConvRows[i]); + idx.push_back(localConvRows[i].event_idx); + } + localOffsets.push_back(static_cast(localStream.tellp())); + std::string localStr = localStream.str(); + + MPI_Comm conv_comm = static_cast(world); + if (world.rank() == 0) { + std::vector > all; + for (size_t i = 0; i < idx.size(); i++) { + std::string row = localStr.substr(localOffsets[i], + localOffsets[i+1] - localOffsets[i]); + all.push_back(std::make_pair(idx[i], row)); + } + for (int p = 1; p < world.size(); p++) { + int n = 0; + MPI_Recv(&n, 1, MPI_INT, p, 10, conv_comm, MPI_STATUS_IGNORE); + if (n <= 0) continue; + std::vector remIdx(n), remOff(n + 1); + MPI_Recv(&remIdx[0], n, MPI_INT, p, 11, conv_comm, MPI_STATUS_IGNORE); + MPI_Recv(&remOff[0], n + 1, MPI_INT, p, 12, conv_comm, + MPI_STATUS_IGNORE); + int total = remOff[n]; + std::string buf(total, '\0'); + if (total > 0) { + MPI_Recv(&buf[0], total, MPI_CHAR, p, 13, conv_comm, + MPI_STATUS_IGNORE); + } + for (int i = 0; i < n; i++) { + all.push_back(std::make_pair( + remIdx[i], + buf.substr(remOff[i], remOff[i+1] - remOff[i]))); + } + } + std::sort(all.begin(), all.end()); + std::string convFile = outputFile + "_convergence.csv"; + std::ofstream cout(convFile.c_str(), + std::ios::out | std::ios::trunc); + cout << "event_idx,contingency,type,converged,iterations," + "final_tolerance,max_p_bus,max_p_mismatch,max_q_bus," + "max_q_mismatch,status_code\n"; + for (size_t i = 0; i < all.size(); i++) cout << all[i].second; + cout.close(); + printf("[convergence] wrote %zu rows to %s\n", + all.size(), convFile.c_str()); + } else { + int n = static_cast(idx.size()); + MPI_Send(&n, 1, MPI_INT, 0, 10, conv_comm); + if (n > 0) { + MPI_Send(&idx[0], n, MPI_INT, 0, 11, conv_comm); + MPI_Send(&localOffsets[0], n + 1, MPI_INT, 0, 12, conv_comm); + int total = localOffsets[n]; + if (total > 0) { + MPI_Send(const_cast(localStr.c_str()), total, MPI_CHAR, 0, 13, + conv_comm); + } } - sout.close(); -#endif } } // Print out statistics on contingencies -#ifdef USE_STATBLOCK - int t_stats = timer->createCategory("Write Statistics"); - timer->start(t_stats); - vmag_stats.writeMeanAndRMS("vmag.txt",1,false); - vmag_stats.writeMinAndMax("vmag_mm.txt",1,false); - if (check_Qlim) vmag_stats.writeMaskValueCount("pq_change_cnt.txt",2,false); - vang_stats.writeMeanAndRMS("vang.txt",1,false); - vang_stats.writeMinAndMax("vang_mm.txt",1,false); - pgen_stats.writeMeanAndRMS("pgen.txt",1); - pgen_stats.writeMinAndMax("pgen_mm.txt",1); - qgen_stats.writeMeanAndRMS("qgen.txt",1); - qgen_stats.writeMinAndMax("qgen_mm.txt",1); - pflow_stats.writeMeanAndRMS("pflow.txt",1); - pflow_stats.writeMinAndMax("pflow_mm.txt",1); - pflow_stats.writeMaskValueCount("line_flt_cnt.txt",2); - qflow_stats.writeMeanAndRMS("qflow.txt",1); - qflow_stats.writeMinAndMax("qflow_mm.txt",1); - perf_stats.writeMinAndMax("perf_mm.txt",1); - perf_stats.sumColumnValues("perf_sum.txt",1); - timer->stop(t_stats); -#endif + if (write_stats) { + int t_stats = timer->createCategory("Write Statistics"); + timer->start(t_stats); + vmag_stats->writeMeanAndRMS("vmag.txt",1,false); + vmag_stats->writeMinAndMax("vmag_mm.txt",1,false); + if (check_Qlim) vmag_stats->writeMaskValueCount("pq_change_cnt.txt",2,false); + vang_stats->writeMeanAndRMS("vang.txt",1,false); + vang_stats->writeMinAndMax("vang_mm.txt",1,false); + pgen_stats->writeMeanAndRMS("pgen.txt",1); + pgen_stats->writeMinAndMax("pgen_mm.txt",1); + qgen_stats->writeMeanAndRMS("qgen.txt",1); + qgen_stats->writeMinAndMax("qgen_mm.txt",1); + pflow_stats->writeMeanAndRMS("pflow.txt",1); + pflow_stats->writeMinAndMax("pflow_mm.txt",1); + pflow_stats->writeMaskValueCount("line_flt_cnt.txt",2); + qflow_stats->writeMeanAndRMS("qflow.txt",1); + qflow_stats->writeMinAndMax("qflow_mm.txt",1); + perf_stats->writeMinAndMax("perf_mm.txt",1); + perf_stats->sumColumnValues("perf_sum.txt",1); + timer->stop(t_stats); + } timer->stop(t_total); // If all processors executed at least one task, then print out timing // statistics (this printout does not work if some processors do not define diff --git a/src/applications/data_sets/input/ca/input.euro.xml b/src/applications/data_sets/input/ca/input.euro.xml index a9812f721..a41c70bd1 100644 --- a/src/applications/data_sets/input/ca/input.euro.xml +++ b/src/applications/data_sets/input/ca/input.euro.xml @@ -3,7 +3,6 @@ true contingencies_euro.xml - 1 1.1 0.9 false diff --git a/src/applications/data_sets/input/ca/input.polish.xml b/src/applications/data_sets/input/ca/input.polish.xml index 36c1ef786..c867bcad4 100644 --- a/src/applications/data_sets/input/ca/input.polish.xml +++ b/src/applications/data_sets/input/ca/input.polish.xml @@ -3,7 +3,6 @@ true contingencies_polish.xml - 1 1.1 0.9 false diff --git a/src/applications/data_sets/input/ca/input_118.xml b/src/applications/data_sets/input/ca/input_118.xml index 33ec77bee..805078565 100644 --- a/src/applications/data_sets/input/ca/input_118.xml +++ b/src/applications/data_sets/input/ca/input_118.xml @@ -8,7 +8,6 @@ - 1 1.1 0.9 false diff --git a/src/applications/data_sets/input/ca/input_14.xml b/src/applications/data_sets/input/ca/input_14.xml index 92e816e86..c347607ef 100644 --- a/src/applications/data_sets/input/ca/input_14.xml +++ b/src/applications/data_sets/input/ca/input_14.xml @@ -8,7 +8,6 @@ true true - 1 1.1 0.9 true diff --git a/src/applications/data_sets/input/ca/input_14_auto_n1.xml b/src/applications/data_sets/input/ca/input_14_auto_n1.xml index 65577b442..a45148d7a 100644 --- a/src/applications/data_sets/input/ca/input_14_auto_n1.xml +++ b/src/applications/data_sets/input/ca/input_14_auto_n1.xml @@ -5,7 +5,6 @@ false true true - 1 1.1 0.9 false diff --git a/src/applications/data_sets/input/ca/input_14_filters_example.xml b/src/applications/data_sets/input/ca/input_14_filters_example.xml new file mode 100644 index 000000000..b1ce7bfdd --- /dev/null +++ b/src/applications/data_sets/input/ca/input_14_filters_example.xml @@ -0,0 +1,76 @@ + + + + + false + true + true + 1.1 + 0.9 + true + false + + + csv_delta + ca_IEEE14 + + + + + + + + + + + + + C + + + + IEEE14_ca.raw + 50 + 1.0e-6 + true + + pre_ + + -ksp_type richardson + -pc_type lu + -pc_factor_mat_solver_type superlu_dist + -ksp_max_it 1 + + + + diff --git a/src/applications/data_sets/input/ca/input_14_qlim.xml b/src/applications/data_sets/input/ca/input_14_qlim.xml index 52d38cdcf..9f3f1157f 100644 --- a/src/applications/data_sets/input/ca/input_14_qlim.xml +++ b/src/applications/data_sets/input/ca/input_14_qlim.xml @@ -4,7 +4,6 @@ false contingencies_14.xml - 1 1.1 0.9 true diff --git a/src/applications/data_sets/input/ca/monitor_branches_14.csv b/src/applications/data_sets/input/ca/monitor_branches_14.csv new file mode 100644 index 000000000..2d1568282 --- /dev/null +++ b/src/applications/data_sets/input/ca/monitor_branches_14.csv @@ -0,0 +1,10 @@ +from_bus,to_bus,ckt +# Sample monitor allowlist for the IEEE 14-bus case. +# Format: from_bus,to_bus,ckt (one branch per row) +# Header row is tolerated; '#' starts a line comment; blank lines are skipped. +# Whitespace and Windows line endings are tolerated. +1,2,1 +1,5,1 +2,3,1 +4,7,1 +4,9,1 diff --git a/src/applications/modules/powerflow/pf_app_module.cpp b/src/applications/modules/powerflow/pf_app_module.cpp index bd33584da..d5fcb0a30 100644 --- a/src/applications/modules/powerflow/pf_app_module.cpp +++ b/src/applications/modules/powerflow/pf_app_module.cpp @@ -1156,9 +1156,13 @@ gridpack::powerflow::PFAppModule::collectResults() bres.qLoss = bres.qFrom + bres.qTo; bres.rateA = branch->getBranchRatingA(cktIds[k]); - if (bres.rateA > 0.0) { + // Under the current contingency rating tier (A/B/C with fallback). + double ratePicked = p_factory->pickBranchRating(i, static_cast(k)); + if (ratePicked <= 0.0) ratePicked = bres.rateA; + bres.rateSelected = ratePicked; + if (ratePicked > 0.0) { double maxMVA = (bres.mvaFrom > bres.mvaTo) ? bres.mvaFrom : bres.mvaTo; - bres.loadingPercent = maxMVA / bres.rateA * 100.0; + bres.loadingPercent = maxMVA / ratePicked * 100.0; } else { bres.loadingPercent = 0.0; } @@ -1872,6 +1876,12 @@ void gridpack::powerflow::PFAppModule::useRateB(bool flag) p_factory->useRateB(flag); } +void gridpack::powerflow::PFAppModule::setContingencyRating( + const std::string& rating) +{ + p_factory->setContingencyRating(rating); +} + /** * Suppress all output from power flow module * @param flag if true, suppress printing diff --git a/src/applications/modules/powerflow/pf_app_module.hpp b/src/applications/modules/powerflow/pf_app_module.hpp index 1582dbd92..5d7ab249a 100644 --- a/src/applications/modules/powerflow/pf_app_module.hpp +++ b/src/applications/modules/powerflow/pf_app_module.hpp @@ -436,6 +436,12 @@ class PFAppModule */ void useRateB(bool flag); + /** + * Select rating tier ("A" | "B" | "C") for CA violation checks and + * loadingPercent. A->B->C fallback when the picked tier is zero. + */ + void setContingencyRating(const std::string& rating); + /** * Suppress all output from power flow module * @param flag if true, suppress printing diff --git a/src/applications/modules/powerflow/pf_factory_module.cpp b/src/applications/modules/powerflow/pf_factory_module.cpp index d2597ec61..858e33c8a 100644 --- a/src/applications/modules/powerflow/pf_factory_module.cpp +++ b/src/applications/modules/powerflow/pf_factory_module.cpp @@ -58,6 +58,7 @@ PFFactoryModule::PFFactoryModule(PFFactoryModule::NetworkPtr network) { p_network = network; p_rateB = false; + p_contingencyRating = "A"; p_islandCount = 0; p_hasLoneBus = false; p_qlim_deadband = 0.1; @@ -171,11 +172,15 @@ bool gridpack::powerflow::PFFactoryModule::checkLoneBus(std::ofstream *stream) bool bus_ok = true; char buf[128]; p_saveIsolatedStatus.clear(); + p_loneBusIndices.clear(); for (i=0; igetActiveBus(i)) continue; gridpack::powerflow::PFBus *bus = dynamic_cast (p_network->getBus(i).get()); + // Skip already-isolated buses (e.g. PSS/E type-4) so they are not + // re-flagged as lone on every call. + if (bus->isIsolated()) continue; std::vector > branches; bus->getNeighborBranches(branches); int size = branches.size(); @@ -200,6 +205,7 @@ bool gridpack::powerflow::PFFactoryModule::checkLoneBus(std::ofstream *stream) if (!ok) { sprintf(buf,"\nLone bus %d found\n",bus->getOriginalIndex()); p_saveIsolatedStatus.push_back(bus->isIsolated()); + p_loneBusIndices.push_back(i); bus->setIsolated(true); printf("%s",buf); if (stream != NULL) *stream << buf; @@ -217,42 +223,17 @@ bool gridpack::powerflow::PFFactoryModule::checkLoneBus(std::ofstream *stream) void gridpack::powerflow::PFFactoryModule::clearLoneBus() { p_hasLoneBus = false; - if (p_saveIsolatedStatus.size() == 0) return; - int numBus = p_network->numBuses(); - int i, j, k; - int ncount = 0; - for (i=0; igetActiveBus(i)) continue; + // Restore status of buses marked by the last checkLoneBus call. + for (size_t k = 0; k < p_loneBusIndices.size(); k++) { + int i = p_loneBusIndices[k]; gridpack::powerflow::PFBus *bus = dynamic_cast (p_network->getBus(i).get()); - std::vector > branches; - bus->getNeighborBranches(branches); - int size = branches.size(); - bool ok = true; - if (size == 0) { - ok = false; - } - if (ok) { - ok = false; - for (j=0; j status = - dynamic_cast - (branches[j].get())->getLineStatus(); - int nlines = status.size(); - for (k=0; kgetOriginalIndex()); - bus->setIsolated(p_saveIsolatedStatus[ncount]); - ncount++; - } + printf("\nLone bus %d reset\n", bus->getOriginalIndex()); + bus->setIsolated(p_saveIsolatedStatus[k]); } + p_loneBusIndices.clear(); + p_saveIsolatedStatus.clear(); } /** @@ -721,53 +702,38 @@ bool gridpack::powerflow::PFFactoryModule::checkLineOverloadViolations() int numBranch = p_network->numBranches(); int i; bool branch_ok = true; + // p_rateB is legacy rtpr; treat it as "rating tier B" for one call. + std::string savedRating = p_contingencyRating; + if (p_rateB) p_contingencyRating = "B"; for (i=0; igetActiveBranch(i)) { gridpack::powerflow::PFBranch *branch = dynamic_cast (p_network->getBranch(i).get()); - // Loop over all lines in the branch and choose the smallest rating value int nlines; p_network->getBranchData(i)->getValue(BRANCH_NUM_ELEMENTS,&nlines); std::vector tags = branch->getLineTags(); - double rate; for (int k = 0; kgetIgnore(tags[k])) { - bool foundRating=false; - if (p_rateB) { - if (p_network->getBranchData(i)->getValue(BRANCH_RATING_B,&rate,k)) { - foundRating = true; - } else { - if (p_network->getBranchData(i)->getValue(BRANCH_RATING_A,&rate,k)) { - foundRating = true; - } - } - } else { - if (p_network->getBranchData(i)->getValue(BRANCH_RATING_A,&rate,k)) { - foundRating = true; - } - } - if (foundRating) { - if (rate > 0.0) { - gridpack::ComplexType s = branch->getComplexPower(tags[k]); - double pq = abs(s); - if (pq > rate) { - branch_ok = false; - gridpack::powerflow::PFFactoryModule::Violation violation; - violation.bus_violation = false; - violation.line_violation = true; - violation.bus1 = branch->getBus1OriginalIndex(); - violation.bus2 = branch->getBus2OriginalIndex(); - strncpy(violation.tag,tags[k].c_str(),2); - violation.tag[2] = '\0'; - p_violations.push_back(violation); - } - } - } + if (branch->getIgnore(tags[k])) continue; + double rate = pickBranchRating(i, k); + if (rate <= 0.0) continue; + gridpack::ComplexType s = branch->getComplexPower(tags[k]); + double pq = abs(s); + if (pq > rate) { + branch_ok = false; + gridpack::powerflow::PFFactoryModule::Violation violation; + violation.bus_violation = false; + violation.line_violation = true; + violation.bus1 = branch->getBus1OriginalIndex(); + violation.bus2 = branch->getBus2OriginalIndex(); + strncpy(violation.tag,tags[k].c_str(),2); + violation.tag[2] = '\0'; + p_violations.push_back(violation); } } } } + p_contingencyRating = savedRating; return checkTrue(branch_ok); } @@ -782,52 +748,34 @@ bool gridpack::powerflow::PFFactoryModule::checkLineOverloadViolations(int area) int numBranch = p_network->numBranches(); int i; bool branch_ok = true; + std::string savedRating = p_contingencyRating; + if (p_rateB) p_contingencyRating = "B"; for (i=0; igetActiveBranch(i)) { gridpack::powerflow::PFBranch *branch = dynamic_cast (p_network->getBranch(i).get()); - // get buses at either end gridpack::powerflow::PFBus *bus1 = dynamic_cast (branch->getBus1().get()); gridpack::powerflow::PFBus *bus2 = dynamic_cast (branch->getBus2().get()); - // Loop over all lines in the branch and choose the smallest rating value - if (bus1->getArea() == area || bus2->getArea() == area) { - int nlines; - p_network->getBranchData(i)->getValue(BRANCH_NUM_ELEMENTS,&nlines); - std::vector tags = branch->getLineTags(); - double rate; - for (int k = 0; kgetIgnore(tags[k])) { - bool foundRating=false; - if (p_rateB) { - if (p_network->getBranchData(i)->getValue(BRANCH_RATING_B,&rate,k)) { - foundRating = true; - } else { - if (p_network->getBranchData(i)->getValue(BRANCH_RATING_A,&rate,k)) { - foundRating = true; - } - } - } else { - if (p_network->getBranchData(i)->getValue(BRANCH_RATING_A,&rate,k)) { - foundRating = true; - } - } - if (foundRating) { - if (rate > 0.0) { - gridpack::ComplexType s = branch->getComplexPower(tags[k]); - double pq = abs(s); - if (pq > rate) branch_ok = false; - } - } - } - } + if (bus1->getArea() != area && bus2->getArea() != area) continue; + int nlines; + p_network->getBranchData(i)->getValue(BRANCH_NUM_ELEMENTS,&nlines); + std::vector tags = branch->getLineTags(); + for (int k = 0; kgetIgnore(tags[k])) continue; + double rate = pickBranchRating(i, k); + if (rate <= 0.0) continue; + gridpack::ComplexType s = branch->getComplexPower(tags[k]); + double pq = abs(s); + if (pq > rate) branch_ok = false; } } } + p_contingencyRating = savedRating; return checkTrue(branch_ok); } @@ -1504,6 +1452,41 @@ void gridpack::powerflow::PFFactoryModule::useRateB(bool flag) } } +/** + * Select rating tier for overload checks. + */ +void gridpack::powerflow::PFFactoryModule::setContingencyRating( + const std::string& rating) +{ + if (rating == "A" || rating == "B" || rating == "C") { + p_contingencyRating = rating; + } else { + p_contingencyRating = "A"; + } + // Keep PFBranch's serialWrite("flow",...) denominator in sync so the .out + // file's loading% matches _violations.csv / JSON loading_percent. + gridpack::powerflow::PFBranch::setContingencyRating(p_contingencyRating); +} + +/** + * Rating for one line element under the current contingency tier + * with A->B->C fallback when the picked tier is zero/missing. + */ +double gridpack::powerflow::PFFactoryModule::pickBranchRating( + int branchLocalIdx, int elemIdx) const +{ + double a = 0.0, b = 0.0, c = 0.0; + p_network->getBranchData(branchLocalIdx)->getValue(BRANCH_RATING_A, &a, elemIdx); + p_network->getBranchData(branchLocalIdx)->getValue(BRANCH_RATING_B, &b, elemIdx); + p_network->getBranchData(branchLocalIdx)->getValue(BRANCH_RATING_C, &c, elemIdx); + if (p_contingencyRating == "A") return a; + if (p_contingencyRating == "B") return (b > 0.0) ? b : a; + // "C" + if (c > 0.0) return c; + if (b > 0.0) return b; + return a; +} + /** * Check switched shunt violations and adjust shunt B values. * For buses with SWREM != 0, resolves remote bus voltage via getLocalBusIndices. diff --git a/src/applications/modules/powerflow/pf_factory_module.hpp b/src/applications/modules/powerflow/pf_factory_module.hpp index 880855175..20fddae61 100644 --- a/src/applications/modules/powerflow/pf_factory_module.hpp +++ b/src/applications/modules/powerflow/pf_factory_module.hpp @@ -365,10 +365,23 @@ class PFFactoryModule */ void useRateB(bool flag); + /** + * Select which rating tier drives overload checks. 'A', 'B', 'C'; + * A->B->C fallback when the picked tier is zero/missing. + */ + void setContingencyRating(const std::string& rating); + std::string getContingencyRating() const { return p_contingencyRating; } + + /** + * Rating for one branch element under the current tier + fallback. + */ + double pickBranchRating(int branchLocalIdx, int elemIdx) const; + private: NetworkPtr p_network; std::vector p_saveIsolatedStatus; + std::vector p_loneBusIndices; std::vector p_saveIslandIsolatedStatus; // For island detection std::vector p_islandIsolatedBusIndices; // Local indices of buses isolated due to islanding int p_islandCount; // Number of islands detected @@ -380,6 +393,7 @@ class PFFactoryModule std::vector p_violations; bool p_rateB; + std::string p_contingencyRating; // "A" | "B" | "C" (default "A") double p_qlim_deadband; // Q deadband (Mvar) for PV->PQ switch }; diff --git a/src/applications/powerflow/pf_main.cpp b/src/applications/powerflow/pf_main.cpp index 28b9f4c8e..f3a08d39c 100644 --- a/src/applications/powerflow/pf_main.cpp +++ b/src/applications/powerflow/pf_main.cpp @@ -76,6 +76,7 @@ int main(int argc, char **argv) // Parse qlim flag (default: true - enforce reactive power limits) bool qlim = cursor->get("qlim", true); + double qlim_deadband = cursor->get("qlimDeadband", 0.1); // Set flags BEFORE creating network // This must be called before readNetwork() for it to take effect @@ -85,6 +86,7 @@ int main(int argc, char **argv) gridpack::powerflow::PFBus::setInitStartMode(gridpack::powerflow::INIT_START_WARM); } gridpack::powerflow::PFBus::setQlim(qlim); + gridpack::powerflow::PFBus::setQlimDeadband(qlim_deadband); // setup and run powerflow calculation boost::shared_ptr diff --git a/src/parallel/communicator.cpp b/src/parallel/communicator.cpp index bd3d2ae98..25aa17753 100644 --- a/src/parallel/communicator.cpp +++ b/src/parallel/communicator.cpp @@ -186,6 +186,13 @@ Communicator::divide(int nsize) const { int nprocs(size()); int me(rank()); + // Fast path: when the result would contain every rank in this comm, return + // a copy that shares the existing GA process group handle. Skips a redundant + // GA_Pgroup_create (whose handle dispatches through a slower path than GA's + // built-in world group on some MPI stacks). + if (nsize >= nprocs) { + return *this; + } // find out how many communicators need to be created int ngrp = nprocs/nsize; if (ngrp*nsize < nprocs) ngrp++; diff --git a/src/parser/PTI33_parser.hpp b/src/parser/PTI33_parser.hpp index 59116e6ab..2315fc2f4 100644 --- a/src/parser/PTI33_parser.hpp +++ b/src/parser/PTI33_parser.hpp @@ -273,13 +273,13 @@ class PTI33_parser : public BasePTIParser<_network> multi_section_parser.parse(p_istream,p_branchData); gridpack::parser::ZoneParser33 zone_parser(&p_busMap, &p_nameMap, &p_branchMap); - zone_parser.parse(p_istream); + zone_parser.parse(p_istream,p_network_data); gridpack::parser::InterAreaParser33 interarea_parser(&p_busMap, &p_nameMap, &p_branchMap); interarea_parser.parse(p_istream); gridpack::parser::OwnerParser33 owner_parser(&p_busMap, &p_nameMap, &p_branchMap); - owner_parser.parse(p_istream); + owner_parser.parse(p_istream,p_network_data); gridpack::parser::FACTSParser33 facts_parser(&p_busMap, &p_nameMap, &p_branchMap); facts_parser.parse(p_istream); diff --git a/src/parser/PTI34_parser.hpp b/src/parser/PTI34_parser.hpp index 248218f52..5d918a021 100644 --- a/src/parser/PTI34_parser.hpp +++ b/src/parser/PTI34_parser.hpp @@ -286,13 +286,13 @@ class PTI34_parser : public BasePTIParser<_network> multi_section_parser.parse(p_istream,p_branchData); gridpack::parser::ZoneParser33 zone_parser(&p_busMap, &p_nameMap, &p_branchMap); - zone_parser.parse(p_istream); + zone_parser.parse(p_istream,p_network_data); gridpack::parser::InterAreaParser33 interarea_parser(&p_busMap, &p_nameMap, &p_branchMap); interarea_parser.parse(p_istream); gridpack::parser::OwnerParser33 owner_parser(&p_busMap, &p_nameMap, &p_branchMap); - owner_parser.parse(p_istream); + owner_parser.parse(p_istream,p_network_data); gridpack::parser::FACTSParser33 facts_parser(&p_busMap, &p_nameMap, &p_branchMap); facts_parser.parse(p_istream); diff --git a/src/parser/PTI35_parser.hpp b/src/parser/PTI35_parser.hpp index d14dbc898..b2368297b 100644 --- a/src/parser/PTI35_parser.hpp +++ b/src/parser/PTI35_parser.hpp @@ -281,13 +281,13 @@ class PTI35_parser : public BasePTIParser<_network> multi_section_parser.parse(p_istream,p_branchData); gridpack::parser::ZoneParser33 zone_parser(&p_busMap, &p_nameMap, &p_branchMap); - zone_parser.parse(p_istream); + zone_parser.parse(p_istream,p_network_data); gridpack::parser::InterAreaParser33 interarea_parser(&p_busMap, &p_nameMap, &p_branchMap); interarea_parser.parse(p_istream); gridpack::parser::OwnerParser33 owner_parser(&p_busMap, &p_nameMap, &p_branchMap); - owner_parser.parse(p_istream); + owner_parser.parse(p_istream,p_network_data); gridpack::parser::FACTSParser33 facts_parser(&p_busMap, &p_nameMap, &p_branchMap); facts_parser.parse(p_istream); diff --git a/src/parser/PTI36_parser.hpp b/src/parser/PTI36_parser.hpp index 3b5e7d47c..d90735cfe 100644 --- a/src/parser/PTI36_parser.hpp +++ b/src/parser/PTI36_parser.hpp @@ -289,13 +289,13 @@ class PTI36_parser : public BasePTIParser<_network> multi_section_parser.parse(p_istream,p_branchData); gridpack::parser::ZoneParser33 zone_parser(&p_busMap, &p_nameMap, &p_branchMap); - zone_parser.parse(p_istream); + zone_parser.parse(p_istream,p_network_data); gridpack::parser::InterAreaParser33 interarea_parser(&p_busMap, &p_nameMap, &p_branchMap); interarea_parser.parse(p_istream); gridpack::parser::OwnerParser33 owner_parser(&p_busMap, &p_nameMap, &p_branchMap); - owner_parser.parse(p_istream); + owner_parser.parse(p_istream,p_network_data); gridpack::parser::FACTSParser33 facts_parser(&p_busMap, &p_nameMap, &p_branchMap); facts_parser.parse(p_istream); diff --git a/src/parser/block_parsers/owner_parser33.cpp b/src/parser/block_parsers/owner_parser33.cpp index e7ed5eb76..807b1823a 100644 --- a/src/parser/block_parsers/owner_parser33.cpp +++ b/src/parser/block_parsers/owner_parser33.cpp @@ -34,34 +34,37 @@ gridpack::parser::OwnerParser33::~OwnerParser33(void) } /** - * parse owner block. Currently does not store data + * parse owner block * @param stream input stream that feeds lines from RAW file + * @param p_network_data data collection object to store parameters from RAW file */ void gridpack::parser::OwnerParser33::parse( - gridpack::stream::InputStream &stream) + gridpack::stream::InputStream &stream, + boost::shared_ptr &p_network_data) { std::string line; stream.nextLine(line); //this should be the first line of the block + int ncnt = 0; while(test_end(line)) { -#if 0 std::vector split_line; if (check_comment(line)) { stream.nextLine(line); continue; } + this->cleanComment(line); split_line = this->splitPSSELine(line); - std::vector owner_instance; - gridpack::component::DataCollection data; - data.addValue(OWNER_NUMBER, atoi(split_line[0].c_str())); - owner_instance.push_back(data); + if (split_line.size() >= 2) { + // OWNER_NUMBER "I" integer + p_network_data->addValue(OWNER_NUMBER, atoi(split_line[0].c_str()), ncnt); - data.addValue(OWNER_NAME, split_line[1].c_str()); - owner_instance.push_back(data); + // OWNER_NAME "OWNAM" string + p_network_data->addValue(OWNER_NAME, split_line[1].c_str(), ncnt); + ncnt++; + } - owner.push_back(owner_instance); -#endif stream.nextLine(line); } + p_network_data->addValue(OWNER_TOTAL, ncnt); } diff --git a/src/parser/block_parsers/owner_parser33.hpp b/src/parser/block_parsers/owner_parser33.hpp index 05b7f19fb..7c510f9e3 100644 --- a/src/parser/block_parsers/owner_parser33.hpp +++ b/src/parser/block_parsers/owner_parser33.hpp @@ -35,11 +35,13 @@ class OwnerParser33 : public BaseBlockParser { virtual ~OwnerParser33(void); /** - * parse owner block. Currently does not store data + * parse owner block * @param stream input stream that feeds lines from RAW file + * @param p_network_data data collection object to store parameters */ void parse( - gridpack::stream::InputStream &stream); + gridpack::stream::InputStream &stream, + boost::shared_ptr &p_network_data); }; } // parser diff --git a/src/parser/block_parsers/transformer_parser34.cpp b/src/parser/block_parsers/transformer_parser34.cpp index 32d08a336..af42576e8 100644 --- a/src/parser/block_parsers/transformer_parser34.cpp +++ b/src/parser/block_parsers/transformer_parser34.cpp @@ -633,14 +633,20 @@ void gridpack::parser::TransformerParser34::parse( /* * type: float - * BRANCH_RATE1-12 + * BRANCH_RATE1-12; also mirror the first three into + * BRANCH_RATING_A/B/C so downstream consumers (PFBranch::load, + * overload checks, JSON export) see the ratings that other PTI + * dialects populate. */ - p_branchData[l_idx]->addValue(BRANCH_RATE1, - atof(split_line3[3].c_str()),nelems); - p_branchData[l_idx]->addValue(BRANCH_RATE2, - atof(split_line3[4].c_str()),nelems); - p_branchData[l_idx]->addValue(BRANCH_RATE3, - atof(split_line3[5].c_str()),nelems); + double rate1 = atof(split_line3[3].c_str()); + double rate2 = atof(split_line3[4].c_str()); + double rate3 = atof(split_line3[5].c_str()); + p_branchData[l_idx]->addValue(BRANCH_RATE1, rate1, nelems); + p_branchData[l_idx]->addValue(BRANCH_RATE2, rate2, nelems); + p_branchData[l_idx]->addValue(BRANCH_RATE3, rate3, nelems); + p_branchData[l_idx]->addValue(BRANCH_RATING_A, rate1, nelems); + p_branchData[l_idx]->addValue(BRANCH_RATING_B, rate2, nelems); + p_branchData[l_idx]->addValue(BRANCH_RATING_C, rate3, nelems); p_branchData[l_idx]->addValue(BRANCH_RATE4, atof(split_line3[6].c_str()),nelems); p_branchData[l_idx]->addValue(BRANCH_RATE5, diff --git a/src/parser/block_parsers/transformer_parser35.cpp b/src/parser/block_parsers/transformer_parser35.cpp index 514d69f81..2bdbc4407 100644 --- a/src/parser/block_parsers/transformer_parser35.cpp +++ b/src/parser/block_parsers/transformer_parser35.cpp @@ -634,14 +634,18 @@ void gridpack::parser::TransformerParser35::parse( /* * type: float - * BRANCH_RATE1-12 - */ - p_branchData[l_idx]->addValue(BRANCH_RATE1, - atof(split_line3[3].c_str()),nelems); - p_branchData[l_idx]->addValue(BRANCH_RATE2, - atof(split_line3[4].c_str()),nelems); - p_branchData[l_idx]->addValue(BRANCH_RATE3, - atof(split_line3[5].c_str()),nelems); + * BRANCH_RATE1-12; also mirror the first three into + * BRANCH_RATING_A/B/C for downstream consumers. + */ + double rate1 = atof(split_line3[3].c_str()); + double rate2 = atof(split_line3[4].c_str()); + double rate3 = atof(split_line3[5].c_str()); + p_branchData[l_idx]->addValue(BRANCH_RATE1, rate1, nelems); + p_branchData[l_idx]->addValue(BRANCH_RATE2, rate2, nelems); + p_branchData[l_idx]->addValue(BRANCH_RATE3, rate3, nelems); + p_branchData[l_idx]->addValue(BRANCH_RATING_A, rate1, nelems); + p_branchData[l_idx]->addValue(BRANCH_RATING_B, rate2, nelems); + p_branchData[l_idx]->addValue(BRANCH_RATING_C, rate3, nelems); p_branchData[l_idx]->addValue(BRANCH_RATE4, atof(split_line3[6].c_str()),nelems); p_branchData[l_idx]->addValue(BRANCH_RATE5, diff --git a/src/parser/block_parsers/zone_parser33.cpp b/src/parser/block_parsers/zone_parser33.cpp index 20cccaad6..a8f511eb1 100644 --- a/src/parser/block_parsers/zone_parser33.cpp +++ b/src/parser/block_parsers/zone_parser33.cpp @@ -34,19 +34,38 @@ gridpack::parser::ZoneParser33::~ZoneParser33(void) } /** - * parse zone block. Currently does not store data + * parse zone block * @param stream input stream that feeds lines from RAW file + * @param p_network_data data collection object to store parameters from RAW file */ void gridpack::parser::ZoneParser33::parse( - gridpack::stream::InputStream &stream) + gridpack::stream::InputStream &stream, + boost::shared_ptr &p_network_data) { std::string line; stream.nextLine(line); //this should be the first line of the block + int ncnt = 0; while(test_end(line)) { - // TODO: parse something here + std::vector split_line; + if (check_comment(line)) { + stream.nextLine(line); + continue; + } + this->cleanComment(line); + split_line = this->splitPSSELine(line); + + if (split_line.size() >= 2) { + // ZONE_NUMBER "I" integer + p_network_data->addValue(ZONE_NUMBER, atoi(split_line[0].c_str()), ncnt); + + // ZONE_NAME "ZONAME" string + p_network_data->addValue(ZONE_NAME, split_line[1].c_str(), ncnt); + ncnt++; + } + stream.nextLine(line); } - + p_network_data->addValue(ZONE_TOTAL, ncnt); } diff --git a/src/parser/block_parsers/zone_parser33.hpp b/src/parser/block_parsers/zone_parser33.hpp index 3d6e836b8..25ec5c793 100644 --- a/src/parser/block_parsers/zone_parser33.hpp +++ b/src/parser/block_parsers/zone_parser33.hpp @@ -35,11 +35,13 @@ class ZoneParser33 : public BaseBlockParser { virtual ~ZoneParser33(void); /** - * parse zone block. Currently does not store data + * parse zone block * @param stream input stream that feeds lines from RAW file + * @param p_network_data data collection object to store parameters */ void parse( - gridpack::stream::InputStream &stream); + gridpack::stream::InputStream &stream, + boost::shared_ptr &p_network_data); }; } // parser diff --git a/src/parser/variable_defs/misc_defs.hpp b/src/parser/variable_defs/misc_defs.hpp index 146021c6e..fa547838f 100644 --- a/src/parser/variable_defs/misc_defs.hpp +++ b/src/parser/variable_defs/misc_defs.hpp @@ -152,15 +152,23 @@ // ZONE DATA +/** + * Total number of zone fields + * type: integer + */ +#define ZONE_TOTAL "ZONE_TOTAL" + /** * Zone Number * type: integer + * indexed */ #define ZONE_NUMBER "ZONE_NUMBER" /** * Zone Name * type: string + * indexed */ #define ZONE_NAME "ZONE_NAME" @@ -193,15 +201,23 @@ // OWNER +/** + * Total number of owner fields + * type: integer + */ +#define OWNER_TOTAL "OWNER_TOTAL" + /** * Owner number * type: integer + * indexed */ #define OWNER_NUMBER "OWNER_NUMBER" /** * Owner name - * type: integer + * type: string + * indexed */ #define OWNER_NAME "OWNER_NAME" diff --git a/src/utilities/results_exporter.cpp b/src/utilities/results_exporter.cpp index b55a0f790..b8e241826 100644 --- a/src/utilities/results_exporter.cpp +++ b/src/utilities/results_exporter.cpp @@ -144,6 +144,7 @@ void ResultsExporter::writeBranchesJSON(std::ostream& out, out << indent << " \"mva_from\": " << std::setprecision(4) << br.mvaFrom << ",\n"; out << indent << " \"mva_to\": " << std::setprecision(4) << br.mvaTo << ",\n"; out << indent << " \"rate_a_mva\": " << std::setprecision(4) << br.rateA << ",\n"; + out << indent << " \"rate_selected_mva\": " << std::setprecision(4) << br.rateSelected << ",\n"; out << indent << " \"loading_percent\": " << std::setprecision(2) << br.loadingPercent << "\n"; out << indent << " }"; if (i + 1 < branches.size()) { @@ -183,6 +184,52 @@ void ResultsExporter::writeGeneratorsJSON(std::ostream& out, out << indent << "]"; } +// ------------------------------------------------------------- +// writeViolationsJSON +// ------------------------------------------------------------- +void ResultsExporter::writeViolationsJSON(std::ostream& out, + const std::vector& brViols, + const std::vector& vViols, + const std::string& indent) +{ + out << std::fixed; + out << indent << "\"violations\": {\n"; + out << indent << " \"branches\": ["; + for (size_t i = 0; i < brViols.size(); ++i) { + const BranchViolation& v = brViols[i]; + out << (i == 0 ? "\n" : ",\n"); + out << indent << " {" + << "\"from_bus\": " << v.fromBus + << ", \"to_bus\": " << v.toBus + << ", \"circuit_id\": \"" << escapeJSON(v.circuitId) << "\"" + << ", \"mva\": " << std::setprecision(4) << v.mva + << ", \"rate_mva\": " << std::setprecision(4) << v.rate + << ", \"loading_percent\": " << std::setprecision(2) << v.loadingPercent + << ", \"base_mva\": " << std::setprecision(4) << v.baseMva + << ", \"delta_mva\": " << std::setprecision(4) << v.deltaMva + << ", \"severity\": \"" << escapeJSON(v.severity) << "\"" + << "}"; + } + if (!brViols.empty()) out << "\n" << indent << " "; + out << "],\n"; + out << indent << " \"voltages\": ["; + for (size_t i = 0; i < vViols.size(); ++i) { + const VoltageViolation& v = vViols[i]; + out << (i == 0 ? "\n" : ",\n"); + out << indent << " {" + << "\"bus_id\": " << v.busId + << ", \"v_pu\": " << std::setprecision(6) << v.vPu + << ", \"limit_low\": " << std::setprecision(4) << v.limitLow + << ", \"limit_high\": " << std::setprecision(4) << v.limitHigh + << ", \"deviation_pu\": " << std::setprecision(6) << v.deviationPu + << ", \"severity\": \"" << escapeJSON(v.severity) << "\"" + << "}"; + } + if (!vViols.empty()) out << "\n" << indent << " "; + out << "]\n"; + out << indent << "}"; +} + // ------------------------------------------------------------- // writePFJSON // ------------------------------------------------------------- @@ -221,23 +268,7 @@ void ResultsExporter::writeCAJSON(std::ofstream& out, out << " \"contingencies\": [\n"; for (size_t i = 0; i < r.contingencies.size(); ++i) { - const ContingencyResult& ct = r.contingencies[i]; - out << " {\n"; - out << " \"name\": \"" << escapeJSON(ct.name) << "\",\n"; - out << " \"type\": \"" << escapeJSON(ct.type) << "\",\n"; - out << " \"has_voltage_violation\": " << (ct.hasVoltageViolation ? "true" : "false") << ",\n"; - out << " \"has_branch_violation\": " << (ct.hasBranchViolation ? "true" : "false") << ",\n"; - out << " \"solution\": {\n"; - writeConvergenceJSON(out, ct.solution.convergence, " "); - out << ",\n"; - writeBusesJSON(out, ct.solution.buses, " "); - out << ",\n"; - writeBranchesJSON(out, ct.solution.branches, " "); - out << ",\n"; - writeGeneratorsJSON(out, ct.solution.generators, " "); - out << "\n"; - out << " }\n"; - out << " }"; + writeContingencyResultJSON(out, r.contingencies[i]); if (i + 1 < r.contingencies.size()) { out << ","; } @@ -320,7 +351,7 @@ void ResultsExporter::writePFCSV(const std::string& basename, out << "from_bus,to_bus,circuit_id," << "p_from_mw,q_from_mvar,p_to_mw,q_to_mvar," << "p_loss_mw,q_loss_mvar," - << "mva_from,mva_to,rate_a_mva,loading_percent\n"; + << "mva_from,mva_to,rate_a_mva,rate_selected_mva,loading_percent\n"; } for (size_t i = 0; i < r.branches.size(); ++i) { @@ -340,6 +371,7 @@ void ResultsExporter::writePFCSV(const std::string& basename, << std::setprecision(4) << br.mvaFrom << "," << std::setprecision(4) << br.mvaTo << "," << std::setprecision(4) << br.rateA << "," + << std::setprecision(4) << br.rateSelected << "," << std::setprecision(2) << br.loadingPercent << "\n"; } out.close(); @@ -485,6 +517,8 @@ void ResultsExporter::writeContingencyResultJSON(std::ostream& out, << (ct.hasVoltageViolation ? "true" : "false") << ",\n"; out << " \"has_branch_violation\": " << (ct.hasBranchViolation ? "true" : "false") << ",\n"; + writeViolationsJSON(out, ct.branchViolations, ct.voltageViolations, " "); + out << ",\n"; out << " \"solution\": {\n"; writeConvergenceJSON(out, ct.solution.convergence, " "); out << ",\n"; diff --git a/src/utilities/results_exporter.hpp b/src/utilities/results_exporter.hpp index 6ad70bb91..49a940fd3 100644 --- a/src/utilities/results_exporter.hpp +++ b/src/utilities/results_exporter.hpp @@ -58,8 +58,9 @@ struct BranchResult { double qLoss; // MVAr double mvaFrom; // MVA double mvaTo; // MVA - double rateA; // MVA - double loadingPercent; // max(|S_from|,|S_to|)/rateA * 100 + double rateA; // MVA (always rate-A) + double rateSelected; // MVA under the picked contingency rating tier + double loadingPercent; // max(|S_from|,|S_to|)/rateSelected * 100 }; struct GeneratorResult { @@ -95,12 +96,35 @@ struct PowerFlowResults { std::vector generators; }; +struct BranchViolation { + int fromBus; + int toBus; + std::string circuitId; + double mva; // max(|S_from|,|S_to|) + double rate; // MVA under selected tier + double loadingPercent; // mva / rate * 100 + double baseMva; // base-case max MVA on same element (0 if unknown) + double deltaMva; // mva - baseMva + std::string severity; // "warning" (<105%), "critical" (>=105%) +}; + +struct VoltageViolation { + int busId; + double vPu; + double limitLow; + double limitHigh; + double deviationPu; // signed: v - nearest limit + std::string severity; +}; + struct ContingencyResult { std::string name; std::string type; // "branch" or "generator" PowerFlowResults solution; bool hasVoltageViolation; bool hasBranchViolation; + std::vector branchViolations; + std::vector voltageViolations; }; struct ContingencyAnalysisResults { @@ -192,6 +216,11 @@ class ResultsExporter { const std::vector& gens, const std::string& indent); + static void writeViolationsJSON(std::ostream& out, + const std::vector& brViols, + const std::vector& vViols, + const std::string& indent); + static void writePFJSON(std::ofstream& out, const PowerFlowResults& r); static void writeCAJSON(std::ofstream& out,