From fe3f8eb09c248ce95641bb35fda6cf787e153090 Mon Sep 17 00:00:00 2001 From: "yousu.chen@pnnl.gov" Date: Thu, 19 Feb 2026 18:34:09 -0800 Subject: [PATCH] fix: update transformer parser and pf_components for ZIP load model --- .../components/pf_matrix/pf_components.cpp | 132 ++++++++++++++++-- .../components/pf_matrix/pf_components.hpp | 15 ++ .../block_parsers/transformer_parser33.cpp | 13 +- .../block_parsers/transformer_parser34.cpp | 13 +- .../block_parsers/transformer_parser35.cpp | 13 +- 5 files changed, 162 insertions(+), 24 deletions(-) diff --git a/src/applications/components/pf_matrix/pf_components.cpp b/src/applications/components/pf_matrix/pf_components.cpp index a2b9f64d..ec6b7b04 100755 --- a/src/applications/components/pf_matrix/pf_components.cpp +++ b/src/applications/components/pf_matrix/pf_components.cpp @@ -28,6 +28,10 @@ * capability (Qmax) to share reactive power among generators * @date 2026-02-02 * + * @updated Yousu Chen + * - Implemented voltage-dependent ZIP load model + * @date 2026-02-19 + * * @brief Methods used in power flow application * * @@ -301,13 +305,17 @@ bool gridpack::powerflow::PFBus::chkQlim(void) return false; } - // Calculate total load Q + // Calculate total load Q (constant-power part) double ql = 0.0; for (int i = 0; i < p_lstatus.size(); i++) { if (p_lstatus[i] == 1) { ql += p_ql[i]; } } + // Add voltage-dependent Q load (IQ*V - YQ*V^2) + double pzip, qzip; + getZIPLoadPower(p_v, pzip, qzip); + ql += qzip; // Required Q from generators = Q_injection + Q_load double Q_required = p_Qinj * p_sbase + ql; @@ -626,6 +634,14 @@ void gridpack::powerflow::PFBus::load( p_pb.clear(); p_pl.clear(); p_ql.clear(); + p_ip.clear(); + p_iq.clear(); + p_yp.clear(); + p_yq.clear(); + p_saveIp.clear(); + p_saveIq.clear(); + p_saveYp.clear(); + p_saveYq.clear(); p_lstatus.clear(); p_lid.clear(); @@ -815,14 +831,19 @@ void gridpack::powerflow::PFBus::load( p_load = p_load && data->getValue(LOAD_YQ, &yq,i); p_load = p_load && data->getValue(LOAD_STATUS, &lstatus,i); if (p_load) { - /* Combine constant P, constant I, constant Y loads */ - pl = pl + ip + yp; - ql = ql - iq - yq; - + /* Store constant-power load only; ZIP components stored separately */ p_pl.push_back(pl); p_savePl.push_back(pl); p_ql.push_back(ql); p_saveQl.push_back(ql); + p_ip.push_back(ip); + p_iq.push_back(iq); + p_yp.push_back(yp); + p_yq.push_back(yq); + p_saveIp.push_back(ip); + p_saveIq.push_back(iq); + p_saveYp.push_back(yp); + p_saveYq.push_back(yq); p_lstatus.push_back(lstatus); std::string id("-1"); data->getValue(LOAD_ID,&id,i); @@ -1044,6 +1065,10 @@ double gridpack::powerflow::PFBus::getTotalGenOutput() pl += p_pl[i]; } } + // Add voltage-dependent P load (IP*V + YP*V^2) + double pzip, qzip; + getZIPLoadPower(p_v, pzip, qzip); + pl += pzip; totalPgen = p_Pinj * p_sbase + pl; } else { // For non-slack buses, sum scheduled generator outputs @@ -1308,6 +1333,13 @@ bool gridpack::powerflow::PFBus::serialWrite(char *string, const int bufsize, if (p_lstatus[i]) pl += p_pl[i]; if (p_lstatus[i]) ql += p_ql[i]; } + // Add voltage-dependent ZIP load at solved voltage + { + double pzip, qzip; + getZIPLoadPower(p_v, pzip, qzip); + pl += pzip; + ql += qzip; + } sprintf(sbuf,"%8d, %4d, %16.8f, %16.8f,",getOriginalIndex(),p_type, pl/p_sbase,ql/p_sbase); int len = strlen(sbuf); @@ -1416,6 +1448,13 @@ bool gridpack::powerflow::PFBus::serialWrite(char *string, const int bufsize, ql += p_ql[i]; } } + // Add voltage-dependent ZIP load at solved voltage + { + double pzip, qzip; + getZIPLoadPower(p_v, pzip, qzip); + pl += pzip; + ql += qzip; + } for (i=0; isetValue("GENERATOR_PF_PGEN",rval,i)) { @@ -1913,17 +1966,29 @@ int gridpack::powerflow::PFBus::getZone() int gridpack::powerflow::PFBus::diagonalJacobianValues(double *rvals) { if (!isIsolated()) { + // Compute ZIP Jacobian contributions: d(pzip)/dV and d(qzip)/dV + double dpzip_dV = 0.0, dqzip_dV = 0.0; + for (int i = 0; i < p_lstatus.size(); i++) { + if (p_lstatus[i] == 1) { + dpzip_dV += p_ip[i] + 2.0 * p_yp[i] * p_v; + dqzip_dV += p_iq[i] - 2.0 * p_yq[i] * p_v; + } + } #ifdef LARGE_MATRIX if (!getReferenceBus()) { - rvals[0] = -p_Qinj - p_ybusi * p_v *p_v; - rvals[1] = p_Pinj - p_ybusr * p_v *p_v; - rvals[2] = p_Pinj / p_v + p_ybusr * p_v; - rvals[3] = p_Qinj / p_v - p_ybusi * p_v; + rvals[0] = -p_Qinj - p_ybusi * p_v *p_v; + rvals[1] = p_Pinj - p_ybusr * p_v *p_v; + rvals[2] = p_Pinj / p_v + p_ybusr * p_v; + rvals[3] = p_Qinj / p_v - p_ybusi * p_v; // Fix up matrix elements if bus is PV bus if (p_isPV) { rvals[1] = 0.0; rvals[2] = 0.0; rvals[3] = 1.0; + } else { + // Add ZIP load derivatives for PQ buses + rvals[2] += dpzip_dV / p_sbase; + rvals[3] += dqzip_dV / p_sbase; } return 4; } else { @@ -1935,14 +2000,16 @@ int gridpack::powerflow::PFBus::diagonalJacobianValues(double *rvals) } #else if (!getReferenceBus() && !p_isPV) { - rvals[0] = -p_Qinj - p_ybusi * p_v *p_v; - rvals[1] = p_Pinj - p_ybusr * p_v *p_v; - rvals[2] = p_Pinj / p_v + p_ybusr * p_v; - rvals[3] = p_Qinj / p_v - p_ybusi * p_v; - // Fix up matrix elements if bus is PV bus + rvals[0] = -p_Qinj - p_ybusi * p_v *p_v; + rvals[1] = p_Pinj - p_ybusr * p_v *p_v; + rvals[2] = p_Pinj / p_v + p_ybusr * p_v; + rvals[3] = p_Qinj / p_v - p_ybusi * p_v; + // Add ZIP load derivatives for PQ buses + rvals[2] += dpzip_dV / p_sbase; + rvals[3] += dqzip_dV / p_sbase; return 4; } else if (!getReferenceBus() && p_isPV) { - rvals[0] = -p_Qinj - p_ybusi * p_v *p_v; + rvals[0] = -p_Qinj - p_ybusi * p_v *p_v; return 1; } else { return 0; @@ -1994,6 +2061,11 @@ int gridpack::powerflow::PFBus::rhsValues(double *rvals) p_Qinj = Q; P -= p_P0; Q -= p_Q0; + // Add voltage-dependent ZIP load contribution + double pzip, qzip; + getZIPLoadPower(p_v, pzip, qzip); + P += pzip / p_sbase; + Q += qzip / p_sbase; rvals[0] = P; #ifdef LARGE_MATRIX if (!p_isPV) { @@ -2185,6 +2257,10 @@ void gridpack::powerflow::PFBus::scaleLoadPower(std::string tag, double value) if (p_lid[i] == tag && p_lstatus[i] == 1) { p_pl[i] = value*p_pl[i]; p_ql[i] = value*p_ql[i]; + p_ip[i] = value*p_ip[i]; + p_iq[i] = value*p_iq[i]; + p_yp[i] = value*p_yp[i]; + p_yq[i] = value*p_yq[i]; break; } } @@ -2203,6 +2279,10 @@ void gridpack::powerflow::PFBus::resetPower() for (i=0; i p_pl, p_ql,p_ip,p_iq,p_yp,p_yq; std::vector p_savePl; std::vector p_saveQl; + std::vector p_saveIp, p_saveIq, p_saveYp, p_saveYq; std::vector p_lstatus; std::vector p_lid; double p_sbase; @@ -598,6 +612,7 @@ class PFBus & p_pt & p_pb & p_pl & p_ql & p_ip & p_iq & p_yp & p_yq & p_savePl & p_saveQl + & p_saveIp & p_saveIq & p_saveYp & p_saveYq & p_lstatus & p_lid & p_sbase & p_Pinj & p_Qinj diff --git a/src/parser/block_parsers/transformer_parser33.cpp b/src/parser/block_parsers/transformer_parser33.cpp index 38cb425c..611d03b9 100644 --- a/src/parser/block_parsers/transformer_parser33.cpp +++ b/src/parser/block_parsers/transformer_parser33.cpp @@ -9,6 +9,8 @@ * Author: Bruce Palmer * Updated 3-winding transformer: Nov. 17th, 2025 * Author: Yousu Chen + * Fixed CZ parameter handling for impedance conversion: Feb. 19th, 2026 + * Author: Yousu Chen */ #include "transformer_parser33.hpp" @@ -431,8 +433,9 @@ void gridpack::parser::TransformerParser33::parse( * type: integer * TRANSFORMER_CZ */ + int cz = atoi(split_line[5].c_str()); p_branchData[l_idx]->addValue(TRANSFORMER_CZ, - atoi(split_line[5].c_str()),nelems); + cz,nelems); /* * type: integer @@ -568,7 +571,9 @@ void gridpack::parser::TransformerParser33::parse( double rval = atof(split_line2[0].c_str()); p_branchData[l_idx]->addValue(TRANSFORMER_R1_2,rval,nelems); rval = rval * windv2 * windv2; // need to consider the wnd2 ratio to the req of the transformer - if (sbase2 == p_case_sbase || sbase2 == 0.0) { + // CZ=1: impedance on system base (no conversion needed) + // CZ=2: impedance on transformer SBASE (convert to system base) + if (cz == 1 || sbase2 == p_case_sbase || sbase2 == 0.0) { p_branchData[l_idx]->addValue(BRANCH_R,rval,nelems); } else { rval = rval*p_case_sbase/sbase2; @@ -583,7 +588,9 @@ void gridpack::parser::TransformerParser33::parse( rval = atof(split_line2[1].c_str()); p_branchData[l_idx]->addValue(TRANSFORMER_X1_2,rval,nelems); rval = rval * windv2 * windv2; // need to consider the wnd2 ratio to the xeq of the transformer - if (sbase2 == p_case_sbase || sbase2 == 0.0) { + // CZ=1: impedance on system base (no conversion needed) + // CZ=2: impedance on transformer SBASE (convert to system base) + if (cz == 1 || sbase2 == p_case_sbase || sbase2 == 0.0) { p_branchData[l_idx]->addValue(BRANCH_X,rval,nelems); } else { rval = rval*p_case_sbase/sbase2; diff --git a/src/parser/block_parsers/transformer_parser34.cpp b/src/parser/block_parsers/transformer_parser34.cpp index 2cf1ad09..1da142e0 100644 --- a/src/parser/block_parsers/transformer_parser34.cpp +++ b/src/parser/block_parsers/transformer_parser34.cpp @@ -9,6 +9,8 @@ * Author: Bruce Palmer * Updated 3-winding transformer: Nov. 17th, 2025 * Author: Yousu Chen + * Fixed CZ parameter handling for impedance conversion: Feb. 19th, 2026 + * Author: Yousu Chen */ #include "transformer_parser34.hpp" #include "ga.h" @@ -468,8 +470,9 @@ void gridpack::parser::TransformerParser34::parse( * type: integer * TRANSFORMER_CZ */ + int cz = atoi(split_line[5].c_str()); p_branchData[l_idx]->addValue(TRANSFORMER_CZ, - atoi(split_line[5].c_str()),nelems); + cz,nelems); /* * type: integer @@ -603,7 +606,9 @@ void gridpack::parser::TransformerParser34::parse( double rval = atof(split_line2[0].c_str()); p_branchData[l_idx]->addValue(TRANSFORMER_R1_2,rval,nelems); rval = rval * windv2 * windv2; // need to consider the wnd2 ratio to the req of the transformer - if (sbase2 == p_case_sbase || sbase2 == 0.0) { + // CZ=1: impedance on system base (no conversion needed) + // CZ=2: impedance on transformer SBASE (convert to system base) + if (cz == 1 || sbase2 == p_case_sbase || sbase2 == 0.0) { p_branchData[l_idx]->addValue(BRANCH_R,rval,nelems); } else { rval = rval*p_case_sbase/sbase2; @@ -618,7 +623,9 @@ void gridpack::parser::TransformerParser34::parse( rval = atof(split_line2[1].c_str()); p_branchData[l_idx]->addValue(TRANSFORMER_X1_2,rval,nelems); rval = rval * windv2 * windv2; // need to consider the wnd2 ratio to the xeq of the transformer - if (sbase2 == p_case_sbase || sbase2 == 0.0) { + // CZ=1: impedance on system base (no conversion needed) + // CZ=2: impedance on transformer SBASE (convert to system base) + if (cz == 1 || sbase2 == p_case_sbase || sbase2 == 0.0) { p_branchData[l_idx]->addValue(BRANCH_X,rval,nelems); } else { rval = rval*p_case_sbase/sbase2; diff --git a/src/parser/block_parsers/transformer_parser35.cpp b/src/parser/block_parsers/transformer_parser35.cpp index 71b19389..c60b908e 100644 --- a/src/parser/block_parsers/transformer_parser35.cpp +++ b/src/parser/block_parsers/transformer_parser35.cpp @@ -9,6 +9,8 @@ * Author: Bruce Palmer * Updated 3-winding transformer: Nov. 17th, 2025 * Author: Yousu Chen + * Fixed CZ parameter handling for impedance conversion: Feb. 19th, 2026 + * Author: Yousu Chen */ #include "transformer_parser35.hpp" @@ -467,8 +469,9 @@ void gridpack::parser::TransformerParser35::parse( * type: integer * TRANSFORMER_CZ */ + int cz = atoi(split_line[5].c_str()); p_branchData[l_idx]->addValue(TRANSFORMER_CZ, - atoi(split_line[5].c_str()),nelems); + cz,nelems); /* * type: integer @@ -603,7 +606,9 @@ void gridpack::parser::TransformerParser35::parse( p_branchData[l_idx]->addValue(TRANSFORMER_R1_2,rval,nelems); rval = rval * windv2 * windv2; // need to consider the wnd2 ratio // to the req of the transformer - if (sbase2 == p_case_sbase || sbase2 == 0.0) { + // CZ=1: impedance on system base (no conversion needed) + // CZ=2: impedance on transformer SBASE (convert to system base) + if (cz == 1 || sbase2 == p_case_sbase || sbase2 == 0.0) { p_branchData[l_idx]->addValue(BRANCH_R,rval,nelems); } else { rval = rval*p_case_sbase/sbase2; @@ -619,7 +624,9 @@ void gridpack::parser::TransformerParser35::parse( p_branchData[l_idx]->addValue(TRANSFORMER_X1_2,rval,nelems); rval = rval * windv2 * windv2; // need to consider the wnd2 ratio // to the xeq of the transformer - if (sbase2 == p_case_sbase || sbase2 == 0.0) { + // CZ=1: impedance on system base (no conversion needed) + // CZ=2: impedance on transformer SBASE (convert to system base) + if (cz == 1 || sbase2 == p_case_sbase || sbase2 == 0.0) { p_branchData[l_idx]->addValue(BRANCH_X,rval,nelems); } else { rval = rval*p_case_sbase/sbase2;