Skip to content

fix: update transformer parser and pf_components for ZIP load model #248

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
chenyousu wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

fix: update transformer parser and pf_components for ZIP load model #248

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
132 changes: 117 additions & 15 deletions src/applications/components/pf_matrix/pf_components.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
*
*
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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();

Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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; i<ngen; i++) {
double pval = 0.0;
double qval = 0.0;
Expand Down Expand Up @@ -1623,6 +1662,13 @@ void gridpack::powerflow::PFBus::saveData(
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; i<ngen; i++) {
if(getReferenceBus()) {
rval = p_pFac[i]*(p_Pinj+pl/p_sbase);
Expand Down Expand Up @@ -1745,6 +1791,13 @@ void gridpack::powerflow::PFBus::saveDataAlsotoOrg(
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; i<ngen; i++) {
rval = p_pFac[i]*(p_Pinj+pl/p_sbase);
if (!data->setValue("GENERATOR_PF_PGEN",rval,i)) {
Expand Down Expand Up @@ -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 {
Expand All @@ -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;
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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;
}
}
Expand All @@ -2203,6 +2279,10 @@ void gridpack::powerflow::PFBus::resetPower()
for (i=0; i<p_nload; i++) {
p_pl[i] = p_savePl[i];
p_ql[i] = p_saveQl[i];
p_ip[i] = p_saveIp[i];
p_iq[i] = p_saveIq[i];
p_yp[i] = p_saveYp[i];
p_yq[i] = p_saveYq[i];
}
}

Expand Down Expand Up @@ -2285,6 +2365,28 @@ void gridpack::powerflow::PFBus::setScale(double scale)
p_rtpr_scale = scale;
}

/**
* Compute total voltage-dependent ZIP load power at given voltage
* PSS/E convention:
* P_load = PL + IP*V + YP*V^2 (all terms add to real load)
* Q_load = QL + IQ*V - YQ*V^2 (IQ=inductive adds, YQ=capacitive subtracts)
* This returns only the voltage-dependent parts:
* pzip = sum(IP*V + YP*V^2)
* qzip = sum(IQ*V - YQ*V^2)
*/
void gridpack::powerflow::PFBus::getZIPLoadPower(double V,
double &pzip, double &qzip) const
{
pzip = 0.0;
qzip = 0.0;
for (int i = 0; i < p_lstatus.size(); i++) {
if (p_lstatus[i] == 1) {
pzip += p_ip[i] * V + p_yp[i] * V * V;
qzip += p_iq[i] * V - p_yq[i] * V * V;
}
}
}

/**
* Simple constructor
*/
Expand Down
15 changes: 15 additions & 0 deletions src/applications/components/pf_matrix/pf_components.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,10 @@
* capability (Qmax) to share reactive power among generators
* @date 2026-02-02
*
* @updated Yousu Chen
* - Added getZIPLoadPower() declaration for voltage-dependent ZIP load model
* @date 2026-02-19
*
* @brief
*
*
Expand Down Expand Up @@ -488,6 +492,15 @@ class PFBus
*/
void setScale(double scale);

/**
* Compute total voltage-dependent ZIP load power at given voltage
* P_zip = sum(IP*V + YP*V^2), Q_zip = sum(IQ*V + YQ*V^2) for active loads
* @param V voltage magnitude
* @param pzip (output) total P contribution from constant-I and constant-Y loads
* @param qzip (output) total Q contribution from constant-I and constant-Y loads
*/
void getZIPLoadPower(double V, double &pzip, double &qzip) const;

/**
* Set the initial start mode for power flow solver
* @param mode INIT_START_WARM (default): use voltage values from raw file
Expand Down Expand Up @@ -542,6 +555,7 @@ class PFBus
std::vector<double> p_pl, p_ql,p_ip,p_iq,p_yp,p_yq;
std::vector<double> p_savePl;
std::vector<double> p_saveQl;
std::vector<double> p_saveIp, p_saveIq, p_saveYp, p_saveYq;
std::vector<int> p_lstatus;
std::vector<std::string> p_lid;
double p_sbase;
Expand Down Expand Up @@ -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
Expand Down
13 changes: 10 additions & 3 deletions src/parser/block_parsers/transformer_parser33.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand Down
13 changes: 10 additions & 3 deletions src/parser/block_parsers/transformer_parser34.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand Down
Loading