diff --git a/datafiles/case9/case9_dcline.m b/datafiles/case9/case9_dcline.m new file mode 100644 index 000000000..0f4bdb3e1 --- /dev/null +++ b/datafiles/case9/case9_dcline.m @@ -0,0 +1,80 @@ +function mpc = case9_dcline +%T_CASE9_DCLINE Same as T_CASE9_OPFV2 with addition of DC line data. +% Please see CASEFORMAT for details on the case file format. +% +% See also: TOGGLE_DCLINE, IDX_DCLINE. + +% MATPOWER + +%% MATPOWER Case Format : Version 2 +mpc.version = '2'; + +%%----- Power Flow Data -----%% +%% system MVA base +mpc.baseMVA = 100; + +%% bus data +% bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin +mpc.bus = [ + 1 3 0 0 0 0 1 1 0 345 1 1.1 0.9; + 2 2 0 0 0 0 1 1 0 345 1 1.1 0.9; + 30 2 0 0 0 0 1 1 0 345 1 1.1 0.9; + 4 1 0 0 0 0 1 1 0 345 1 1.1 0.9; + 5 1 90 30 0 0 1 1 0 345 1 1.1 0.9; + 6 1 0 0 0 0 1 1 0 345 1 1.1 0.9; + 7 1 100 35 0 0 1 1 0 345 1 1.1 0.9; + 8 1 0 0 0 0 1 1 0 345 1 1.1 0.9; + 9 1 125 50 0 0 1 1 0 345 1 1.1 0.9; +]; + +%% generator data +% bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf +mpc.gen = [ + 1 0 0 300 -300 1 100 1 250 90 0 0 0 0 0 0 0 0 0 0 0; + 2 163 0 300 -300 1 100 1 300 10 0 200 -20 20 -10 10 0 0 0 0 0; + 30 85 0 300 -300 1 100 1 270 10 0 200 -30 30 -15 15 0 0 0 0 0; +]; + +%% branch data +% fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax +mpc.branch = [ + 1 4 0 0.0576 0 0 250 250 0 0 1 -360 2.48; + 4 5 0.017 0.092 0.158 0 250 250 0 0 1 -360 360; + 5 6 0.039 0.17 0.358 150 150 150 0 0 1 -360 360; + 30 6 0 0.0586 0 0 300 300 0 0 1 -360 360; + 6 7 0.0119 0.1008 0.209 40 150 150 0 0 1 -360 360; + 7 8 0.0085 0.072 0.149 250 250 250 0 0 1 -360 360; + 8 2 0 0.0625 0 250 250 250 0 0 1 -360 360; + 8 9 0.032 0.161 0.306 250 250 250 0 0 1 -360 360; + 9 4 0.01 0.085 0.176 250 250 250 0 0 1 -2 360; +]; + +%%----- OPF Data -----%% +%% generator cost data +% 1 startup shutdown n x1 y1 ... xn yn +% 2 startup shutdown n c(n-1) ... c0 +mpc.gencost = [ + 2 1500.00 0.00 3 0.11 5 150; + 2 2000.00 0.00 3 0.085 1.2 600; + 2 3000 0 3 0.1225 1 335; +]; + +%%----- DC Line Data ----- +%% Change 30-4 Vf to 1.0 from 1.01 to match with gen set-point voltage +% fbus tbus status Pf Pt Qf Qt Vf Vt Pmin Pmax QminF QmaxF QminT QmaxT loss0 loss1 +mpc.dcline = [ + 30 4 1 10 8.9 0 0 1 1 1 10 -10 10 -10 10 1 0.01; + 7 9 1 2 1.96 0 0 1 0.98 2 10 0 0 0 0 0 0; + 5 8 1 0 0 0 0 1 1 1 10 -10 10 -10 10 0 0; + 5 9 1 10 9.5 0 0 1 0.98 0 10 -10 10 -10 10 0 0.05; +]; + +%% DC line cost data +% 1 startup shutdown n x1 y1 ... xn yn +% 2 startup shutdown n c(n-1) ... c0 +mpc.dclinecost = [ + 2 0 0 2 0 0 0 0 0 0 0 0 0 0; + 2 0 0 2 0 0 0 0 0 0 0 0 0 0; + 2 0 0 2 0 0 0 0 0 0 0 0 0 0; + 2 0 0 2 7.3 0 0 0 0 0 0 0 0 0; +]; diff --git a/include/constants.h b/include/constants.h index e1a27f3d1..621c955d2 100644 --- a/include/constants.h +++ b/include/constants.h @@ -18,7 +18,8 @@ */ #define NLOAD_AT_BUS_MAX 32 /**< Maximum number of loads allowed at a bus */ #define NPHASE 1 /**< Per-phase analysis */ -#define MAX_KV_LEVELS 20 /**< Maximum KV levels */ +#define MAX_KV_LEVELS \ + 200 /**< Maximum KV (voltage) levels allowed in input file */ /* Fuel types for generators */ #define GENFUEL_COAL 0 /* Coal */ diff --git a/include/private/psimpl.h b/include/private/psimpl.h index fea18eb73..1d61cf638 100644 --- a/include/private/psimpl.h +++ b/include/private/psimpl.h @@ -300,6 +300,24 @@ struct _p_PSLINE { PSBUS connbuses[2]; /**< From and to buses */ + /****** For DC lines **********/ + PetscBool isdcline; /**< Is line a DC line? */ + PetscScalar pmin; /**< lower limit on PF (MW flow at "from" end) */ + PetscScalar pmax; /**< upper limit on PF (MW flow at "from" end) */ + PetscScalar Vf; /**ps); CHKERRQ(ierr); + /* set individual load costs if necessary */ PetscFunctionReturn(0); diff --git a/src/opflow/model/dcopf/dcopf.cpp b/src/opflow/model/dcopf/dcopf.cpp index 827bbadc6..fd08d9fc2 100644 --- a/src/opflow/model/dcopf/dcopf.cpp +++ b/src/opflow/model/dcopf/dcopf.cpp @@ -18,6 +18,7 @@ PetscErrorCode OPFLOWSetVariableBounds_DCOPF(OPFLOW opflow, Vec Xl, Vec Xu) { PetscScalar *xl, *xu; PetscInt i; PSBUS bus; + PSLINE line; PSGEN gen; PetscInt loc; @@ -29,6 +30,18 @@ PetscErrorCode OPFLOWSetVariableBounds_DCOPF(OPFLOW opflow, Vec Xl, Vec Xu) { ierr = VecGetArray(Xu, &xu); CHKERRQ(ierr); + for (i = 0; i < opflow->nlinesmon; i++) { + line = &ps->line[opflow->linesmon[i]]; + + if (line->isdcline) { + loc = line->startxdcloc; + + // Bounds on PF + xl[loc] = line->pmin; + xu[loc] = line->pmax; + } + } + for (i = 0; i < ps->nbus; i++) { PetscInt k; @@ -173,12 +186,10 @@ PetscErrorCode OPFLOWSetConstraintBounds_DCOPF(OPFLOW opflow, Vec Gl, Vec Gu) { } } /* Line flow constraint bounds */ - if (!opflow->ignore_lineflow_constraints) { - for (i = 0; i < ps->nline; i++) { - line = &ps->line[i]; - if (!line->status || line->rateA > 1e5) - continue; + for (i = 0; i < opflow->nlinesmon; i++) { + line = &ps->line[opflow->linesmon[i]]; + if (!line->isdcline) { gloc = line->startineqloc; /* Line flow inequality constraints */ gl[gloc] = gl[gloc + 1] = -(line->rateA / ps->MVAbase); @@ -214,10 +225,11 @@ PetscErrorCode OPFLOWSetInitialGuess_DCOPF(OPFLOW opflow, Vec X, Vec Lambda) { PetscErrorCode ierr; PS ps = opflow->ps; const PetscScalar *xl, *xu; - PetscScalar *x; + PetscScalar *x, *lambda; PetscInt i; PSBUS bus; - PetscInt loc; + PSLINE line; + PetscInt loc, gloc; PetscFunctionBegin; @@ -228,6 +240,8 @@ PetscErrorCode OPFLOWSetInitialGuess_DCOPF(OPFLOW opflow, Vec X, Vec Lambda) { CHKERRQ(ierr); ierr = VecGetArrayRead(opflow->Xu, &xu); CHKERRQ(ierr); + ierr = VecGetArray(Lambda, &lambda); + CHKERRQ(ierr); for (i = 0; i < ps->nbus; i++) { PetscInt k; @@ -299,6 +313,26 @@ PetscErrorCode OPFLOWSetInitialGuess_DCOPF(OPFLOW opflow, Vec X, Vec Lambda) { x[loc] = 0.0; } } + + gloc = bus->starteqloc; + lambda[gloc] = bus->mult_pmis; + } + + for (i = 0; i < opflow->nlinesmon; i++) { + line = &ps->line[opflow->linesmon[i]]; + + if (line->isdcline) { + loc = line->startxdcloc; + + if (opflow->initializationtype == OPFLOWINIT_MIDPOINT || + opflow->initializationtype == OPFLOWINIT_FLATSTART) { + x[loc] = 0.5 * (xl[loc] + xu[loc]); + } else if (opflow->initializationtype == OPFLOWINIT_FROMFILE || + opflow->initializationtype == OPFLOWINIT_ACPF || + opflow->initializationtype == OPFLOWINIT_DCOPF) { + x[loc] = PetscMax(line->pmin, PetscMin(line->pf, line->pmax)); + } + } } ierr = VecRestoreArray(X, &x); @@ -307,6 +341,8 @@ PetscErrorCode OPFLOWSetInitialGuess_DCOPF(OPFLOW opflow, Vec X, Vec Lambda) { CHKERRQ(ierr); ierr = VecRestoreArrayRead(opflow->Xu, &xu); CHKERRQ(ierr); + ierr = VecRestoreArray(Lambda, &lambda); + CHKERRQ(ierr); PetscFunctionReturn(0); } @@ -425,31 +461,45 @@ PetscErrorCode OPFLOWComputeEqualityConstraints_DCOPF(OPFLOW opflow, Vec X, busf = connbuses[0]; bust = connbuses[1]; - xlocf = busf->startxVloc; - xloct = bust->startxVloc; + if (!line->isdcline) { + xlocf = busf->startxVloc; + xloct = bust->startxVloc; - thetaf = x[xlocf]; - thetat = x[xloct]; - thetaft = thetaf - thetat; - thetatf = thetat - thetaf; + thetaf = x[xlocf]; + thetat = x[xloct]; + thetaft = thetaf - thetat; + thetatf = thetat - thetaf; - if (bus == busf) { - Pf = Bdc * thetaft + Pshift; + if (bus == busf) { + Pf = Bdc * thetaft + Pshift; - val[0] = Pf; + val[0] = Pf; - ierr = VecSetValues(Ge, 1, row, val, ADD_VALUES); - CHKERRQ(ierr); + ierr = VecSetValues(Ge, 1, row, val, ADD_VALUES); + CHKERRQ(ierr); - flps += 1.0; - } else { - Pt = Bdc * thetatf - Pshift; + flps += 2.0; + } else { + Pt = Bdc * thetatf - Pshift; - val[0] = Pt; + val[0] = Pt; + ierr = VecSetValues(Ge, 1, row, val, ADD_VALUES); + CHKERRQ(ierr); + flps += 2.0; + } + } else if (line->isdcline) { + Pf = x[line->startxdcloc]; + + if (bus == busf) { + val[0] = Pf; + } else { + Pt = Pf - (line->loss0 + line->loss1 * Pf); + val[0] = -Pf; + flps += 3.0; + } ierr = VecSetValues(Ge, 1, row, val, ADD_VALUES); CHKERRQ(ierr); - flps += 1.0; } } @@ -586,10 +636,9 @@ PetscErrorCode OPFLOWComputeEqualityConstraintJacobian_DCOPF(OPFLOW opflow, for (k = 0; k < nconnlines; k++) { line = connlines[k]; + if (!line->status) continue; - Bdc = line->bdc; - // Pshift = line->pshift; /* Get the connected buses to this line */ ierr = PSLINEGetConnectedBuses(line, &connbuses); @@ -597,41 +646,47 @@ PetscErrorCode OPFLOWComputeEqualityConstraintJacobian_DCOPF(OPFLOW opflow, busf = connbuses[0]; bust = connbuses[1]; - // locf = busf->startxVloc; - // loct = bust->startxVloc; - - locglobf = busf->startxVlocglob; - locglobt = bust->startxVlocglob; + if (!line->isdcline) { + locglobf = busf->startxVlocglob; + locglobt = bust->startxVlocglob; - // thetaf = xarr[locf]; - // thetat = xarr[loct]; - // thetaft = thetaf - thetat; - // thetatf = thetat - thetaf; + Bdc = line->bdc; - if (bus == busf) { - col[0] = locglobf; - col[1] = locglobt; - /* dPf_dthetaf */ - val[0] = Bdc; - /*dPf_dthetat */ - val[1] = -Bdc; + if (bus == busf) { + col[0] = locglobf; + col[1] = locglobt; + /* dPf_dthetaf */ + val[0] = Bdc; + /*dPf_dthetat */ + val[1] = -Bdc; - ierr = MatSetValues(Je, 1, row, 2, col, val, ADD_VALUES); - CHKERRQ(ierr); - } else { - col[0] = locglobt; - col[1] = locglobf; + ierr = MatSetValues(Je, 1, row, 2, col, val, ADD_VALUES); + CHKERRQ(ierr); + } else { + col[0] = locglobt; + col[1] = locglobf; - /* dPt_dthetat */ - val[0] = Bdc; - /* dPt_dthetaf */ - val[1] = -Bdc; + /* dPt_dthetat */ + val[0] = Bdc; + /* dPt_dthetaf */ + val[1] = -Bdc; - ierr = MatSetValues(Je, 1, row, 2, col, val, ADD_VALUES); - CHKERRQ(ierr); + ierr = MatSetValues(Je, 1, row, 2, col, val, ADD_VALUES); + CHKERRQ(ierr); + } + } else if (line->isdcline) { + if (bus == busf) { + col[0] = line->startxdcloc; + val[0] = 1.0; + ierr = MatSetValues(Je, 1, row, 1, col, val, ADD_VALUES); + } else { + col[0] = line->startxdcloc; + val[0] = -(1.0 - line->loss1); + ierr = MatSetValues(Je, 1, row, 1, col, val, ADD_VALUES); + flps += 2.0; + } } } - flps += nconnlines * 2; if (opflow->has_gensetpoint) { ierr = PSBUSGetVariableGlobalLocation(bus, &locglob); @@ -741,36 +796,36 @@ PetscErrorCode OPFLOWComputeInequalityConstraints_DCOPF(OPFLOW opflow, Vec X, } if (!opflow->ignore_lineflow_constraints) { - for (i = 0; i < ps->nline; i++) { - line = &ps->line[i]; - if (!line->status || line->rateA > 1e5) - continue; + for (i = 0; i < opflow->nlinesmon; i++) { + line = &ps->line[opflow->linesmon[i]]; - gloc = line->startineqloc; + if (!line->isdcline) { + gloc = line->startineqloc; - Bdc = line->bdc; - Pshift = line->pshift; + Bdc = line->bdc; + Pshift = line->pshift; - ierr = PSLINEGetConnectedBuses(line, &connbuses); - CHKERRQ(ierr); - busf = connbuses[0]; - bust = connbuses[1]; + ierr = PSLINEGetConnectedBuses(line, &connbuses); + CHKERRQ(ierr); + busf = connbuses[0]; + bust = connbuses[1]; - xlocf = busf->startxVloc; - xloct = bust->startxVloc; + xlocf = busf->startxVloc; + xloct = bust->startxVloc; - thetaf = x[xlocf]; - thetat = x[xloct]; - thetaft = thetaf - thetat; - thetatf = thetat - thetaf; + thetaf = x[xlocf]; + thetat = x[xloct]; + thetaft = thetaf - thetat; + thetatf = thetat - thetaf; - Pf = Bdc * thetaft + Pshift; - Pt = Bdc * thetatf - Pshift; + Pf = Bdc * thetaft + Pshift; + Pt = Bdc * thetatf - Pshift; - g[gloc] = Pf; - g[gloc + 1] = Pt; + g[gloc] = Pf; + g[gloc + 1] = Pt; - flps += 4.0; + flps += 4.0; + } } } @@ -859,54 +914,48 @@ PetscErrorCode OPFLOWComputeInequalityConstraintJacobian_DCOPF(OPFLOW opflow, } if (!opflow->ignore_lineflow_constraints) { - for (i = 0; i < ps->nline; i++) { - line = &ps->line[i]; - if (!line->status || line->rateA > 1e5) - continue; - - gloc = line->startineqloc; + for (i = 0; i < opflow->nlinesmon; i++) { + line = &ps->line[opflow->linesmon[i]]; - Bdc = line->bdc; - // Pshift = line->pshift; + if (!line->isdcline) { + gloc = line->startineqloc; - ierr = PSLINEGetConnectedBuses(line, &connbuses); - CHKERRQ(ierr); - busf = connbuses[0]; - bust = connbuses[1]; + Bdc = line->bdc; - xlocf = busf->startxVloc; - xloct = bust->startxVloc; + ierr = PSLINEGetConnectedBuses(line, &connbuses); + CHKERRQ(ierr); + busf = connbuses[0]; + bust = connbuses[1]; - // thetaf = x[xlocf]; - // thetat = x[xloct]; - // thetaft = thetaf - thetat; - // thetatf = thetat - thetaf; + xlocf = busf->startxVloc; + xloct = bust->startxVloc; - dPf_dthetaf = Bdc; - dPf_dthetat = -Bdc; + dPf_dthetaf = Bdc; + dPf_dthetat = -Bdc; - dPt_dthetat = Bdc; - dPt_dthetaf = -Bdc; + dPt_dthetat = Bdc; + dPt_dthetaf = -Bdc; - row[0] = gloc; - col[0] = xlocf; - col[1] = xloct; + row[0] = gloc; + col[0] = xlocf; + col[1] = xloct; - val[0] = dPf_dthetaf; - val[1] = dPf_dthetat; + val[0] = dPf_dthetaf; + val[1] = dPf_dthetat; - ierr = MatSetValues(Ji, 1, row, 2, col, val, ADD_VALUES); - CHKERRQ(ierr); + ierr = MatSetValues(Ji, 1, row, 2, col, val, ADD_VALUES); + CHKERRQ(ierr); - row[0] = gloc + 1; - col[0] = xloct; - col[1] = xlocf; + row[0] = gloc + 1; + col[0] = xloct; + col[1] = xlocf; - val[0] = dPt_dthetat; - val[1] = dPt_dthetaf; + val[0] = dPt_dthetat; + val[1] = dPt_dthetaf; - ierr = MatSetValues(Ji, 1, row, 2, col, val, ADD_VALUES); - CHKERRQ(ierr); + ierr = MatSetValues(Ji, 1, row, 2, col, val, ADD_VALUES); + CHKERRQ(ierr); + } } flps += ps->nline * 2; } @@ -1112,15 +1161,22 @@ PetscErrorCode OPFLOWModelSetNumVariables_DCOPF(OPFLOW opflow, PSLINE line; PetscErrorCode ierr; PetscBool isghost; + PetscInt monidx; PetscFunctionBegin; *nx = 0; - /* No variables for the branches */ - for (i = 0; i < ps->nline; i++) { - line = &ps->line[i]; - branchnvar[i] = line->nx = 0; - *nx += branchnvar[i]; + /* Variables for the branches */ + for (i = 0; i < opflow->nlinesmon; i++) { + monidx = opflow->linesmon[i]; + line = &ps->line[opflow->linesmon[i]]; + + if (line->isdcline) { + branchnvar[monidx] = line->nx = 1; + } else { + branchnvar[monidx] = line->nx = 0; + } + *nx += branchnvar[monidx]; } /* Variables for the buses */ @@ -1243,12 +1299,14 @@ PetscErrorCode OPFLOWModelSetNumConstraints_DCOPF(OPFLOW opflow, } if (!opflow->ignore_lineflow_constraints) { - for (i = 0; i < ps->nline; i++) { - line = &ps->line[i]; - if (line->status && line->rateA < 1e5) { + for (i = 0; i < opflow->nlinesmon; i++) { + line = &ps->line[opflow->linesmon[i]]; + if (!line->isdcline) { *nconineq += 2; /* Real power flow line constraint (from and to bus) */ line->nconineq = 2; line->nconeq = 0; + } else if (line->isdcline) { + line->nconeq = line->nconineq = 0; } } } @@ -1555,35 +1613,49 @@ PetscErrorCode OPFLOWSolutionToPS_DCOPF(OPFLOW opflow) { continue; } - Bdc = line->bdc; - Pshift = line->pshift; + if (!line->isdcline) { + Bdc = line->bdc; + Pshift = line->pshift; - ierr = PSLINEGetConnectedBuses(line, &connbuses); - CHKERRQ(ierr); - busf = connbuses[0]; - bust = connbuses[1]; + ierr = PSLINEGetConnectedBuses(line, &connbuses); + CHKERRQ(ierr); + busf = connbuses[0]; + bust = connbuses[1]; - xlocf = busf->startxVloc; - xloct = bust->startxVloc; + xlocf = busf->startxVloc; + xloct = bust->startxVloc; - thetaf = x[xlocf]; - thetat = x[xloct]; - thetaft = thetaf - thetat; - thetatf = thetat - thetaf; + thetaf = x[xlocf]; + thetat = x[xloct]; + thetaft = thetaf - thetat; + thetatf = thetat - thetaf; - Pf = Bdc * thetaft + Pshift; - Pt = Bdc * thetatf - Pshift; + Pf = Bdc * thetaft + Pshift; + Pt = Bdc * thetatf - Pshift; - line->pf = Pf; - line->qf = 0.0; - line->pt = Pt; - line->qt = 0.0; - line->sf = Pf; - line->st = Pt; + line->pf = Pf; + line->qf = 0.0; + line->pt = Pt; + line->qt = 0.0; + line->sf = Pf; + line->st = Pt; + } else if (line->isdcline) { + Pf = x[line->startxdcloc]; + + Pt = Pf - (line->loss0 + line->loss1 * Pf); + + line->pf = Pf; + line->qf = 0.0; + line->pt = Pt; + line->qt = 0.0; + line->sf = Pf; + line->st = Pt; + } + } - if (opflow->ignore_lineflow_constraints || line->rateA > 1e5) { - line->mult_sf = line->mult_st = 0.0; - } else { + for (i = 0; i < opflow->nlinesmon; i++) { + line = &ps->line[opflow->linesmon[i]]; + if (!line->isdcline) { gloc = line->startineqloc; line->mult_sf = lambdai[gloc]; line->mult_st = lambdai[gloc + 1]; @@ -1697,11 +1769,16 @@ PetscErrorCode OPFLOWModelSetUp_DCOPF(OPFLOW opflow) { } } - for (i = 0; i < ps->nline; i++) { - line = &ps->line[i]; - if (line->status && line->rateA < 1e5) { + for (i = 0; i < opflow->nlinesmon; i++) { + line = &ps->line[opflow->linesmon[i]]; + + if (!line->isdcline) { line->startineqloc = ineqloc; ineqloc += line->nconineq; + } else if (line->isdcline) { + ierr = PSLINEGetVariableLocation(line, &loc); + CHKERRQ(ierr); + line->startxdcloc = loc; } } diff --git a/src/opflow/model/power_bal_polar/pbpol.cpp b/src/opflow/model/power_bal_polar/pbpol.cpp index 2f6bfb9a6..d67e16cbd 100644 --- a/src/opflow/model/power_bal_polar/pbpol.cpp +++ b/src/opflow/model/power_bal_polar/pbpol.cpp @@ -1,7 +1,6 @@ #include "pbpol.h" #include "exago_config.h" #include -/* Testing Orestis's permission to commit. */ PetscErrorCode OPFLOWModelDestroy_PBPOL(OPFLOW opflow) { PetscErrorCode ierr; @@ -33,11 +32,28 @@ PetscErrorCode OPFLOWSetVariableBounds_PBPOL(OPFLOW opflow, Vec Xl, Vec Xu) { for (i = 0; i < opflow->nlinesmon; i++) { line = &ps->line[opflow->linesmon[i]]; - if (opflow->allow_lineflow_violation) { - /* Bounds on slack variables */ - loc = line->startxslackloc; - xl[loc] = xl[loc + 1] = 0.0; - xu[loc] = xu[loc + 1] = PETSC_INFINITY; + + if (!line->isdcline) { + if (opflow->allow_lineflow_violation) { + /* Bounds on slack variables */ + loc = line->startxslackloc; + xl[loc] = xl[loc + 1] = 0.0; + xu[loc] = xu[loc + 1] = PETSC_INFINITY; + } + } else { + loc = line->startxdcloc; + + // Bounds on PF + xl[loc] = line->pmin; + xu[loc] = line->pmax; + + // Bounds on QF + xl[loc + 1] = line->qminf; + xu[loc + 1] = line->qmaxf; + + // Bounds on QT + xl[loc + 2] = line->qmint; + xu[loc + 2] = line->qmaxt; } } @@ -97,7 +113,15 @@ PetscErrorCode OPFLOWSetVariableBounds_PBPOL(OPFLOW opflow, Vec Xl, Vec Xu) { loc = gen->startxpowloc; - xl[loc] = gen->pb; /* PGmin */ + /* If generator is renewable then set the lower + bound to 0 so that it can be curtailed if needed + Note: Do this only if we have positive Pmax (pt) + */ + if (gen->isrenewable && gen->pt > gen->pb) { + xl[loc] = 0.0; + } else { + xl[loc] = gen->pb; /* PGmin */ + } xu[loc] = gen->pt; /* PGmax */ xl[loc + 1] = gen->qb; /* QGmin */ xu[loc + 1] = gen->qt; /* QGmax */ @@ -232,11 +256,13 @@ PetscErrorCode OPFLOWSetConstraintBounds_PBPOL(OPFLOW opflow, Vec Gl, Vec Gu) { for (i = 0; i < opflow->nlinesmon; i++) { line = &ps->line[opflow->linesmon[i]]; - gloc = line->startineqloc; - /* Line flow inequality constraints */ - gl[gloc] = gl[gloc + 1] = 0.0; - gu[gloc] = gu[gloc + 1] = - (line->rateA / ps->MVAbase) * (line->rateA / ps->MVAbase); + if (!line->isdcline) { + gloc = line->startineqloc; + /* Line flow inequality constraints */ + gl[gloc] = gl[gloc + 1] = 0.0; + gu[gloc] = gu[gloc + 1] = + (line->rateA / ps->MVAbase) * (line->rateA / ps->MVAbase); + } } ierr = VecRestoreArray(Gl, &gl); @@ -389,14 +415,33 @@ PetscErrorCode OPFLOWSetInitialGuess_PBPOL(OPFLOW opflow, Vec X, Vec Lambda) { for (i = 0; i < opflow->nlinesmon; i++) { line = &ps->line[opflow->linesmon[i]]; - if (opflow->allow_lineflow_violation) { - loc = line->startxslackloc; - /* Initialize slacks for line flow violations */ - x[loc] = x[loc + 1] = 0.0; + if (!line->isdcline) { + if (opflow->allow_lineflow_violation) { + loc = line->startxslackloc; + /* Initialize slacks for line flow violations */ + x[loc] = x[loc + 1] = 0.0; + } + + gloc = line->startineqloc; + lambdai[gloc] = line->mult_sf; + lambdai[gloc + 1] = line->mult_st; + } else if (line->isdcline) { + loc = line->startxdcloc; + + if (opflow->initializationtype == OPFLOWINIT_MIDPOINT || + opflow->initializationtype == OPFLOWINIT_FLATSTART) { + x[loc] = 0.5 * (xl[loc] + xu[loc]); + x[loc + 1] = 0.5 * (xl[loc + 1] + xu[loc + 1]); + x[loc + 2] = 0.5 * (xl[loc + 2] + xu[loc + 2]); + + } else if (opflow->initializationtype == OPFLOWINIT_FROMFILE || + opflow->initializationtype == OPFLOWINIT_ACPF || + opflow->initializationtype == OPFLOWINIT_DCOPF) { + x[loc] = PetscMax(line->pmin, PetscMin(line->pf, line->pmax)); + x[loc + 1] = PetscMax(line->qminf, PetscMin(line->qf, line->qmaxf)); + x[loc + 2] = PetscMax(line->qmint, PetscMin(line->qt, line->qmaxt)); + } } - gloc = line->startineqloc; - lambdai[gloc] = line->mult_sf; - lambdai[gloc + 1] = line->mult_st; } ierr = VecRestoreArray(X, &x); @@ -534,53 +579,70 @@ PetscErrorCode OPFLOWComputeEqualityConstraints_PBPOL(OPFLOW opflow, Vec X, if (!line->status) continue; - Gff = line->yff[0]; - Bff = line->yff[1]; - Gft = line->yft[0]; - Bft = line->yft[1]; - Gtf = line->ytf[0]; - Btf = line->ytf[1]; - Gtt = line->ytt[0]; - Btt = line->ytt[1]; - ierr = PSLINEGetConnectedBuses(line, &connbuses); CHKERRQ(ierr); busf = connbuses[0]; bust = connbuses[1]; - xlocf = busf->startxVloc; - xloct = bust->startxVloc; - - thetaf = x[xlocf]; - Vmf = x[xlocf + 1]; - thetat = x[xloct]; - Vmt = x[xloct + 1]; - thetaft = thetaf - thetat; - thetatf = thetat - thetaf; - - if (bus == busf) { - Pf = Gff * Vmf * Vmf + - Vmf * Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); - Qf = -Bff * Vmf * Vmf + - Vmf * Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); - - val[0] = Pf; - val[1] = Qf; - ierr = VecSetValues(Ge, 2, row, val, ADD_VALUES); - CHKERRQ(ierr); + if (!line->isdcline) { + Gff = line->yff[0]; + Bff = line->yff[1]; + Gft = line->yft[0]; + Bft = line->yft[1]; + Gtf = line->ytf[0]; + Btf = line->ytf[1]; + Gtt = line->ytt[0]; + Btt = line->ytt[1]; + + xlocf = busf->startxVloc; + xloct = bust->startxVloc; + + thetaf = x[xlocf]; + Vmf = x[xlocf + 1]; + thetat = x[xloct]; + Vmt = x[xloct + 1]; + thetaft = thetaf - thetat; + thetatf = thetat - thetaf; + + if (bus == busf) { + Pf = Gff * Vmf * Vmf + + Vmf * Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); + Qf = -Bff * Vmf * Vmf + + Vmf * Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); + + val[0] = Pf; + val[1] = Qf; + ierr = VecSetValues(Ge, 2, row, val, ADD_VALUES); + CHKERRQ(ierr); - flps += 78.0; - } else { - Pt = Gtt * Vmt * Vmt + - Vmt * Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); - Qt = -Btt * Vmt * Vmt + - Vmt * Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); + flps += 78.0; + } else { + Pt = Gtt * Vmt * Vmt + + Vmt * Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); + Qt = -Btt * Vmt * Vmt + + Vmt * Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); - val[0] = Pt; - val[1] = Qt; + val[0] = Pt; + val[1] = Qt; + ierr = VecSetValues(Ge, 2, row, val, ADD_VALUES); + CHKERRQ(ierr); + flps += 78.0; + } + } else if (line->isdcline) { + Pf = x[line->startxdcloc]; + Qf = x[line->startxdcloc + 1]; + Qt = x[line->startxdcloc + 2]; + + if (bus == busf) { + val[0] = Pf; + val[1] = Qf; + } else { + Pt = Pf - (line->loss0 + line->loss1 * Pf); + val[0] = -Pt; + val[1] = -Qt; + } ierr = VecSetValues(Ge, 2, row, val, ADD_VALUES); CHKERRQ(ierr); - flps += 78.0; } } @@ -744,14 +806,6 @@ PetscErrorCode OPFLOWComputeEqualityConstraintJacobian_PBPOL(OPFLOW opflow, line = connlines[k]; if (!line->status) continue; - Gff = line->yff[0]; - Bff = line->yff[1]; - Gft = line->yft[0]; - Bft = line->yft[1]; - Gtf = line->ytf[0]; - Btf = line->ytf[1]; - Gtt = line->ytt[0]; - Btt = line->ytt[1]; /* Get the connected buses to this line */ ierr = PSLINEGetConnectedBuses(line, &connbuses); @@ -759,71 +813,100 @@ PetscErrorCode OPFLOWComputeEqualityConstraintJacobian_PBPOL(OPFLOW opflow, busf = connbuses[0]; bust = connbuses[1]; - locf = busf->startxVloc; - loct = bust->startxVloc; - - locglobf = busf->startxVlocglob; - locglobt = bust->startxVlocglob; - - thetaf = xarr[locf]; - Vmf = xarr[locf + 1]; - thetat = xarr[loct]; - Vmt = xarr[loct + 1]; - thetaft = thetaf - thetat; - thetatf = thetat - thetaf; - - if (bus == busf) { - col[0] = locglobf; - col[1] = locglobf + 1; - col[2] = locglobt; - col[3] = locglobt + 1; - /* dPf_dthetaf */ - val[0] = Vmf * Vmt * (-Gft * sin(thetaft) + Bft * cos(thetaft)); - /*dPf_dVmf */ - val[1] = - 2 * Gff * Vmf + Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); - /*dPf_dthetat */ - val[2] = Vmf * Vmt * (Gft * sin(thetaft) - Bft * cos(thetaft)); - /* dPf_dVmt */ - val[3] = Vmf * (Gft * cos(thetaft) + Bft * sin(thetaft)); - - /* dQf_dthetaf */ - val[4] = Vmf * Vmt * (Bft * sin(thetaft) + Gft * cos(thetaft)); - /* dQf_dVmf */ - val[5] = - -2 * Bff * Vmf + Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); - /* dQf_dthetat */ - val[6] = Vmf * Vmt * (-Bft * sin(thetaft) - Gft * cos(thetaft)); - /* dQf_dVmt */ - val[7] = Vmf * (-Bft * cos(thetaft) + Gft * sin(thetaft)); - ierr = MatSetValues(Je, 2, row, 4, col, val, ADD_VALUES); - CHKERRQ(ierr); - } else { - col[0] = locglobt; - col[1] = locglobt + 1; - col[2] = locglobf; - col[3] = locglobf + 1; - /* dPt_dthetat */ - val[0] = Vmt * Vmf * (-Gtf * sin(thetatf) + Btf * cos(thetatf)); - /* dPt_dVmt */ - val[1] = - 2 * Gtt * Vmt + Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); - /* dPt_dthetaf */ - val[2] = Vmt * Vmf * (Gtf * sin(thetatf) - Btf * cos(thetatf)); - /* dPt_dVmf */ - val[3] = Vmt * (Gtf * cos(thetatf) + Btf * sin(thetatf)); - - /* dQt_dthetat */ - val[4] = Vmt * Vmf * (Btf * sin(thetatf) + Gtf * cos(thetatf)); - /* dQt_dVmt */ - val[5] = - -2 * Btt * Vmt + Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); - /* dQt_dthetaf */ - val[6] = Vmt * Vmf * (-Btf * sin(thetatf) - Gtf * cos(thetatf)); - /* dQt_dVmf */ - val[7] = Vmt * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); - ierr = MatSetValues(Je, 2, row, 4, col, val, ADD_VALUES); - CHKERRQ(ierr); + if (!line->isdcline) { + Gff = line->yff[0]; + Bff = line->yff[1]; + Gft = line->yft[0]; + Bft = line->yft[1]; + Gtf = line->ytf[0]; + Btf = line->ytf[1]; + Gtt = line->ytt[0]; + Btt = line->ytt[1]; + + locf = busf->startxVloc; + loct = bust->startxVloc; + + locglobf = busf->startxVlocglob; + locglobt = bust->startxVlocglob; + + thetaf = xarr[locf]; + Vmf = xarr[locf + 1]; + thetat = xarr[loct]; + Vmt = xarr[loct + 1]; + thetaft = thetaf - thetat; + thetatf = thetat - thetaf; + + if (bus == busf) { + col[0] = locglobf; + col[1] = locglobf + 1; + col[2] = locglobt; + col[3] = locglobt + 1; + /* dPf_dthetaf */ + val[0] = Vmf * Vmt * (-Gft * sin(thetaft) + Bft * cos(thetaft)); + /*dPf_dVmf */ + val[1] = + 2 * Gff * Vmf + Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); + /*dPf_dthetat */ + val[2] = Vmf * Vmt * (Gft * sin(thetaft) - Bft * cos(thetaft)); + /* dPf_dVmt */ + val[3] = Vmf * (Gft * cos(thetaft) + Bft * sin(thetaft)); + + /* dQf_dthetaf */ + val[4] = Vmf * Vmt * (Bft * sin(thetaft) + Gft * cos(thetaft)); + /* dQf_dVmf */ + val[5] = + -2 * Bff * Vmf + Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); + /* dQf_dthetat */ + val[6] = Vmf * Vmt * (-Bft * sin(thetaft) - Gft * cos(thetaft)); + /* dQf_dVmt */ + val[7] = Vmf * (-Bft * cos(thetaft) + Gft * sin(thetaft)); + ierr = MatSetValues(Je, 2, row, 4, col, val, ADD_VALUES); + CHKERRQ(ierr); + } else { + col[0] = locglobt; + col[1] = locglobt + 1; + col[2] = locglobf; + col[3] = locglobf + 1; + /* dPt_dthetat */ + val[0] = Vmt * Vmf * (-Gtf * sin(thetatf) + Btf * cos(thetatf)); + /* dPt_dVmt */ + val[1] = + 2 * Gtt * Vmt + Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); + /* dPt_dthetaf */ + val[2] = Vmt * Vmf * (Gtf * sin(thetatf) - Btf * cos(thetatf)); + /* dPt_dVmf */ + val[3] = Vmt * (Gtf * cos(thetatf) + Btf * sin(thetatf)); + + /* dQt_dthetat */ + val[4] = Vmt * Vmf * (Btf * sin(thetatf) + Gtf * cos(thetatf)); + /* dQt_dVmt */ + val[5] = + -2 * Btt * Vmt + Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); + /* dQt_dthetaf */ + val[6] = Vmt * Vmf * (-Btf * sin(thetatf) - Gtf * cos(thetatf)); + /* dQt_dVmf */ + val[7] = Vmt * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); + ierr = MatSetValues(Je, 2, row, 4, col, val, ADD_VALUES); + CHKERRQ(ierr); + } + } else if (line->isdcline) { + if (bus == busf) { + col[0] = line->startxdcloc; + val[0] = 1.0; + ierr = MatSetValues(Je, 1, row, 1, col, val, ADD_VALUES); + + col[0] = line->startxdcloc + 1; + val[0] = 1.0; + ierr = MatSetValues(Je, 1, row + 1, 1, col, val, ADD_VALUES); + } else { + col[0] = line->startxdcloc; + val[0] = -(1.0 - line->loss1); + ierr = MatSetValues(Je, 1, row, 1, col, val, ADD_VALUES); + + col[0] = line->startxdcloc + 2; + val[0] = -1.0; + ierr = MatSetValues(Je, 1, row + 1, 1, col, val, ADD_VALUES); + } } } flps += nconnlines * @@ -976,6 +1059,9 @@ PetscErrorCode OPFLOWComputeInequalityConstraints_PBPOL(OPFLOW opflow, Vec X, for (i = 0; i < opflow->nlinesmon; i++) { line = &ps->line[opflow->linesmon[i]]; + if (line->isdcline) + continue; + gloc = line->startineqloc; Gff = line->yff[0]; @@ -1029,7 +1115,6 @@ PetscErrorCode OPFLOWComputeInequalityConstraints_PBPOL(OPFLOW opflow, Vec X, g[gloc] -= xsft_slack; g[gloc + 1] -= xstf_slack; } - flps += 160.0; } } @@ -1175,6 +1260,9 @@ PetscErrorCode OPFLOWComputeInequalityConstraintJacobian_PBPOL(OPFLOW opflow, for (i = 0; i < opflow->nlinesmon; i++) { line = &ps->line[opflow->linesmon[i]]; + if (line->isdcline) + continue; + gloc = line->startineqloc; Gff = line->yff[0]; @@ -1274,7 +1362,7 @@ PetscErrorCode OPFLOWComputeInequalityConstraintJacobian_PBPOL(OPFLOW opflow, ierr = MatSetValues(Ji, 1, row, 4, col, val, ADD_VALUES); CHKERRQ(ierr); - if (opflow->allow_lineflow_violation) { + if (!line->isdcline && opflow->allow_lineflow_violation) { loc = line->startxslackloc; row[0] = gloc; col[0] = loc; @@ -1339,6 +1427,8 @@ PetscErrorCode OPFLOWComputeObjective_PBPOL(OPFLOW opflow, Vec X, for (i = 0; i < opflow->nlinesmon; i++) { line = &ps->line[opflow->linesmon[i]]; + if (line->isdcline) + continue; if (opflow->allow_lineflow_violation) { loc = line->startxslackloc; // Slack variables for from and to side @@ -1377,6 +1467,7 @@ PetscErrorCode OPFLOWComputeObjective_PBPOL(OPFLOW opflow, Vec X, if (opflow->objectivetype == MIN_GEN_COST) { loc = gen->startxpowloc; Pg = x[loc] * ps->MVAbase; + *obj += gen->cost_alpha * Pg * Pg + gen->cost_beta * Pg + gen->cost_gamma; flps += 7.0; @@ -1438,13 +1529,13 @@ PetscErrorCode OPFLOWComputeGradient_PBPOL(OPFLOW opflow, Vec X, Vec grad) { for (i = 0; i < opflow->nlinesmon; i++) { line = &ps->line[opflow->linesmon[i]]; + + if (line->isdcline) + continue; if (opflow->allow_lineflow_violation) { loc = line->startxslackloc; - // ADD GRADIENT HERE df[loc] = opflow->lineflowviolation_penalty; df[loc + 1] = opflow->lineflowviolation_penalty; - - flps += 0.0; } } @@ -1537,24 +1628,31 @@ PetscErrorCode OPFLOWModelSetNumVariables_PBPOL(OPFLOW opflow, PSLINE line; PetscErrorCode ierr; PetscBool isghost; - PetscInt idx; + PetscInt monidx; // Index of monitored line PetscFunctionBegin; *nx = 0; - /* No variables for the branches */ + + /* Variables for the lines */ for (i = 0; i < opflow->nlinesmon; i++) { - idx = opflow->linesmon[i]; - line = &ps->line[idx]; - branchnvar[idx] = line->nx = 0; - if (opflow->allow_lineflow_violation) { - /* Two variables for line flow slacks - - From side flow (Sft) and To side flow (Stf) - */ - branchnvar[idx] += 2; - line->nx += 2; + monidx = opflow->linesmon[i]; + line = &ps->line[monidx]; + + if (line->isdcline) { + branchnvar[monidx] = line->nx = 3; + } else { + branchnvar[monidx] = line->nx = 0; + if (opflow->allow_lineflow_violation) { + /* Two variables for line flow slacks + - From side flow (Sft) and To side flow (Stf) + */ + branchnvar[monidx] += 2; + line->nx += 2; + } } - *nx += branchnvar[idx]; + + *nx += branchnvar[monidx]; } /* Variables for the buses */ @@ -1686,9 +1784,13 @@ PetscErrorCode OPFLOWModelSetNumConstraints_PBPOL(OPFLOW opflow, for (i = 0; i < opflow->nlinesmon; i++) { line = &ps->line[opflow->linesmon[i]]; - *nconineq += 2; /* Number of line flow constraints */ - line->nconineq = 2; - line->nconeq = 0; + if (!line->isdcline) { + *nconineq += 2; /* Number of line flow constraints */ + line->nconineq = 2; + line->nconeq = 0; + } else if (line->isdcline) { + line->nconeq = line->nconineq = 0; + } } } @@ -2216,6 +2318,9 @@ PetscErrorCode OPFLOWComputeInequalityConstraintsHessian_PBPOL(OPFLOW opflow, for (i = 0; i < opflow->nlinesmon; i++) { line = &ps->line[opflow->linesmon[i]]; + if (line->isdcline) + continue; + PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt; Gff = line->yff[0]; Bff = line->yff[1]; @@ -2879,40 +2984,47 @@ PetscErrorCode OPFLOWSolutionToPS_PBPOL(OPFLOW opflow) { continue; } - Gff = line->yff[0]; - Bff = line->yff[1]; - Gft = line->yft[0]; - Bft = line->yft[1]; - Gtf = line->ytf[0]; - Btf = line->ytf[1]; - Gtt = line->ytt[0]; - Btt = line->ytt[1]; + if (!line->isdcline) { + Gff = line->yff[0]; + Bff = line->yff[1]; + Gft = line->yft[0]; + Bft = line->yft[1]; + Gtf = line->ytf[0]; + Btf = line->ytf[1]; + Gtt = line->ytt[0]; + Btt = line->ytt[1]; - ierr = PSLINEGetConnectedBuses(line, &connbuses); - CHKERRQ(ierr); - busf = connbuses[0]; - bust = connbuses[1]; + ierr = PSLINEGetConnectedBuses(line, &connbuses); + CHKERRQ(ierr); + busf = connbuses[0]; + bust = connbuses[1]; - xlocf = busf->startxVloc; - xloct = bust->startxVloc; + xlocf = busf->startxVloc; + xloct = bust->startxVloc; - thetaf = x[xlocf]; - Vmf = x[xlocf + 1]; - thetat = x[xloct]; - Vmt = x[xloct + 1]; - thetaft = thetaf - thetat; - thetatf = thetat - thetaf; + thetaf = x[xlocf]; + Vmf = x[xlocf + 1]; + thetat = x[xloct]; + Vmt = x[xloct + 1]; + thetaft = thetaf - thetat; + thetatf = thetat - thetaf; - Pf = - Gff * Vmf * Vmf + Vmf * Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); - Qf = -Bff * Vmf * Vmf + - Vmf * Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); + Pf = Gff * Vmf * Vmf + + Vmf * Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); + Qf = -Bff * Vmf * Vmf + + Vmf * Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); - Pt = - Gtt * Vmt * Vmt + Vmt * Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); - Qt = -Btt * Vmt * Vmt + - Vmt * Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); + Pt = Gtt * Vmt * Vmt + + Vmt * Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); + Qt = -Btt * Vmt * Vmt + + Vmt * Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); + } else if (line->isdcline) { + Pf = x[line->startxdcloc]; + Qf = x[line->startxdcloc + 1]; + Qt = x[line->startxdcloc + 2]; + Pt = Pf - (line->loss0 + line->loss1 * Pf); + } line->pf = Pf; line->qf = Qf; line->pt = Pt; @@ -2923,6 +3035,8 @@ PetscErrorCode OPFLOWSolutionToPS_PBPOL(OPFLOW opflow) { for (i = 0; i < opflow->nlinesmon; i++) { line = &ps->line[opflow->linesmon[i]]; + if (line->isdcline) + continue; gloc = line->startineqloc; line->mult_sf = lambdai[gloc]; line->mult_st = lambdai[gloc + 1]; @@ -3037,16 +3151,18 @@ PetscErrorCode OPFLOWModelSetUp_PBPOL(OPFLOW opflow) { for (i = 0; i < opflow->nlinesmon; i++) { line = &ps->line[opflow->linesmon[i]]; - - ierr = PSLINEGetVariableLocation(line, &loc); - CHKERRQ(ierr); - /* Set starting location for slack variable */ - if (opflow->allow_lineflow_violation) { - line->startxslackloc = loc; + if (!line->isdcline) { + /* Set starting location for slack variable */ + if (opflow->allow_lineflow_violation) { + line->startxslackloc = loc; + } + line->startineqloc = ineqloc; + ineqloc += line->nconineq; + } else if (line->isdcline) { + ierr = PSLINEGetVariableLocation(line, &loc); + CHKERRQ(ierr); + line->startxdcloc = loc; } - - line->startineqloc = ineqloc; - ineqloc += line->nconineq; } PetscFunctionReturn(0); diff --git a/src/pflow/pflow.cpp b/src/pflow/pflow.cpp index a1d58defb..bc34b9c40 100644 --- a/src/pflow/pflow.cpp +++ b/src/pflow/pflow.cpp @@ -178,7 +178,7 @@ PetscErrorCode PFLOWJacobian(SNES snes, Vec X, Mat J, Mat Jpre, void *ctx) { Vec localX; const PetscScalar *xarr; PetscBool ghostbus; - PetscInt i, k; + PetscInt i, j, k; PetscInt flps = 0; PetscFunctionBegin; @@ -233,13 +233,22 @@ PetscErrorCode PFLOWJacobian(SNES snes, Vec X, Mat J, Mat Jpre, void *ctx) { val[0] = 0.0; val[1] = 2 * Vm * bus->gl; val[2] = 0.0; + val[3] = 0.0; flps += 2; if (bus->ide != PV_BUS) { val[3] = -2 * Vm * bus->bl; /* Partial derivative for shunt contribution */ flps += 2; - } else - val[3] = 1.0; /* Partial derivative of voltage magnitude constraint */ + } else { + for (j = 0; j < bus->ngen; j++) { + PSGEN gen; + ierr = PSBUSGetGen(bus, j, &gen); + CHKERRQ(ierr); + if (gen->status) + val[3] += + 1.0; /* Partial derivative of voltage magnitude constraint */ + } + } ierr = MatSetValues(J, 2, row, 2, col, val, ADD_VALUES); CHKERRQ(ierr); } @@ -257,16 +266,6 @@ PetscErrorCode PFLOWJacobian(SNES snes, Vec X, Mat J, Mat Jpre, void *ctx) { line = connlines[k]; if (!line->status) continue; - PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt; - Gff = line->yff[0]; - Bff = line->yff[1]; - Gft = line->yft[0]; - Bft = line->yft[1]; - Gtf = line->ytf[0]; - Btf = line->ytf[1]; - Gtt = line->ytt[0]; - Btt = line->ytt[1]; - const PSBUS *connbuses; PSBUS busf, bust; PetscInt locglobf, locglobt, locf, loct; @@ -296,62 +295,89 @@ PetscErrorCode PFLOWJacobian(SNES snes, Vec X, Mat J, Mat Jpre, void *ctx) { thetatf = thetat - thetaf; flps += 2; - if (bus == busf) { - if (bus->ide != REF_BUS) { - row[0] = locglobf; - col[0] = locglobf; - col[1] = locglobf + 1; - col[2] = locglobt; - col[3] = locglobt + 1; - val[0] = Vmf * Vmt * (-Gft * sin(thetaft) + Bft * cos(thetaft)); - val[1] = - 2 * Gff * Vmf + Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); - val[2] = Vmf * Vmt * (Gft * sin(thetaft) - Bft * cos(thetaft)); - val[3] = Vmf * (Gft * cos(thetaft) + Bft * sin(thetaft)); - flps += 21 + (4 * EXAGO_FLOPS_SINOP) + (4 * EXAGO_FLOPS_COSOP); - ierr = MatSetValues(J, 1, row, 4, col, val, ADD_VALUES); - CHKERRQ(ierr); - - if (bus->ide != PV_BUS) { - row[0] = locglobf + 1; - val[0] = Vmf * Vmt * (Bft * sin(thetaft) + Gft * cos(thetaft)); - val[1] = -2 * Bff * Vmf + - Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); - val[2] = Vmf * Vmt * (-Bft * sin(thetaft) - Gft * cos(thetaft)); - val[3] = Vmf * (-Bft * cos(thetaft) + Gft * sin(thetaft)); + if (!line->isdcline) { + PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt; + Gff = line->yff[0]; + Bff = line->yff[1]; + Gft = line->yft[0]; + Bft = line->yft[1]; + Gtf = line->ytf[0]; + Btf = line->ytf[1]; + Gtt = line->ytt[0]; + Btt = line->ytt[1]; + + if (bus == busf) { + if (bus->ide != REF_BUS) { + row[0] = locglobf; + col[0] = locglobf; + col[1] = locglobf + 1; + col[2] = locglobt; + col[3] = locglobt + 1; + val[0] = Vmf * Vmt * (-Gft * sin(thetaft) + Bft * cos(thetaft)); + val[1] = + 2 * Gff * Vmf + Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); + val[2] = Vmf * Vmt * (Gft * sin(thetaft) - Bft * cos(thetaft)); + val[3] = Vmf * (Gft * cos(thetaft) + Bft * sin(thetaft)); flps += 21 + (4 * EXAGO_FLOPS_SINOP) + (4 * EXAGO_FLOPS_COSOP); ierr = MatSetValues(J, 1, row, 4, col, val, ADD_VALUES); CHKERRQ(ierr); - } - } - } else { - if (bus->ide != REF_BUS) { - row[0] = locglobt; - col[0] = locglobt; - col[1] = locglobt + 1; - col[2] = locglobf; - col[3] = locglobf + 1; - val[0] = Vmt * Vmf * (-Gtf * sin(thetatf) + Btf * cos(thetatf)); - val[1] = - 2 * Gtt * Vmt + Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); - val[2] = Vmt * Vmf * (Gtf * sin(thetatf) - Btf * cos(thetatf)); - val[3] = Vmt * (Gtf * cos(thetatf) + Btf * sin(thetatf)); - flps += 21 + (4 * EXAGO_FLOPS_SINOP) + (4 * EXAGO_FLOPS_COSOP); - ierr = MatSetValues(J, 1, row, 4, col, val, ADD_VALUES); - CHKERRQ(ierr); - if (bus->ide != PV_BUS) { - row[0] = locglobt + 1; - val[0] = Vmt * Vmf * (Btf * sin(thetatf) + Gtf * cos(thetatf)); - val[1] = -2 * Btt * Vmt + - Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); - val[2] = Vmt * Vmf * (-Btf * sin(thetatf) - Gtf * cos(thetatf)); - val[3] = Vmt * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); + if (bus->ide != PV_BUS) { + row[0] = locglobf + 1; + val[0] = Vmf * Vmt * (Bft * sin(thetaft) + Gft * cos(thetaft)); + val[1] = -2 * Bff * Vmf + + Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); + val[2] = Vmf * Vmt * (-Bft * sin(thetaft) - Gft * cos(thetaft)); + val[3] = Vmf * (-Bft * cos(thetaft) + Gft * sin(thetaft)); + flps += 21 + (4 * EXAGO_FLOPS_SINOP) + (4 * EXAGO_FLOPS_COSOP); + ierr = MatSetValues(J, 1, row, 4, col, val, ADD_VALUES); + CHKERRQ(ierr); + } + } + } else { + if (bus->ide != REF_BUS) { + row[0] = locglobt; + col[0] = locglobt; + col[1] = locglobt + 1; + col[2] = locglobf; + col[3] = locglobf + 1; + val[0] = Vmt * Vmf * (-Gtf * sin(thetatf) + Btf * cos(thetatf)); + val[1] = + 2 * Gtt * Vmt + Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); + val[2] = Vmt * Vmf * (Gtf * sin(thetatf) - Btf * cos(thetatf)); + val[3] = Vmt * (Gtf * cos(thetatf) + Btf * sin(thetatf)); flps += 21 + (4 * EXAGO_FLOPS_SINOP) + (4 * EXAGO_FLOPS_COSOP); ierr = MatSetValues(J, 1, row, 4, col, val, ADD_VALUES); CHKERRQ(ierr); + + if (bus->ide != PV_BUS) { + row[0] = locglobt + 1; + val[0] = Vmt * Vmf * (Btf * sin(thetatf) + Gtf * cos(thetatf)); + val[1] = -2 * Btt * Vmt + + Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); + val[2] = Vmt * Vmf * (-Btf * sin(thetatf) - Gtf * cos(thetatf)); + val[3] = Vmt * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); + flps += 21 + (4 * EXAGO_FLOPS_SINOP) + (4 * EXAGO_FLOPS_COSOP); + ierr = MatSetValues(J, 1, row, 4, col, val, ADD_VALUES); + CHKERRQ(ierr); + } } } + } else if (line->isdcline) { + /* DC line */ + if (bus == busf) { + row[0] = locglobf + 1; + col[0] = locglobf + 1; + val[0] = 1.0; + ierr = MatSetValues(J, 1, row, 1, col, val, ADD_VALUES); + CHKERRQ(ierr); + } else { + row[0] = locglobt + 1; + col[0] = locglobt + 1; + val[0] = 1.0; + ierr = MatSetValues(J, 1, row, 1, col, val, ADD_VALUES); + CHKERRQ(ierr); + } } } } @@ -460,7 +486,7 @@ PetscErrorCode PFLOWFunction(SNES snes, Vec X, Vec F, void *ctx) { farr[loc] -= gen->pg; flps += 1; if (bus->ide == PV_BUS) { - farr[loc + 1] = Vm - gen->vs; + farr[loc + 1] += Vm - gen->vs; flps += 1; } else { farr[loc + 1] -= gen->qg; @@ -496,27 +522,18 @@ PetscErrorCode PFLOWFunction(SNES snes, Vec X, Vec F, void *ctx) { line = connlines[k]; if (!line->status) continue; - PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt; - Gff = line->yff[0]; - Bff = line->yff[1]; - Gft = line->yft[0]; - Bft = line->yft[1]; - Gtf = line->ytf[0]; - Btf = line->ytf[1]; - Gtt = line->ytt[0]; - Btt = line->ytt[1]; const PSBUS *connbuses; PSBUS busf, bust; - PetscInt locf, loct; - PetscScalar thetaf, Vmf, thetat, Vmt, thetaft, thetatf; - /* Get the connected buses to this line */ ierr = PSLINEGetConnectedBuses(line, &connbuses); CHKERRQ(ierr); busf = connbuses[0]; bust = connbuses[1]; + PetscInt locf, loct; + PetscScalar thetaf, Vmf, thetat, Vmt, thetaft, thetatf; + ierr = PSBUSGetVariableLocation(busf, &locf); CHKERRQ(ierr); ierr = PSBUSGetVariableLocation(bust, &loct); @@ -525,36 +542,62 @@ PetscErrorCode PFLOWFunction(SNES snes, Vec X, Vec F, void *ctx) { Vmf = xarr[locf + 1]; thetat = xarr[loct]; Vmt = xarr[loct + 1]; - thetaft = thetaf - thetat; - thetatf = thetat - thetaf; - if (bus == busf) { - if (busf->ide != REF_BUS) { - farr[locf] += Gff * Vmf * Vmf + - Vmf * Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); - flps += 9 + EXAGO_FLOPS_SINOP + EXAGO_FLOPS_COSOP; - } - if (busf->ide != PV_BUS && busf->ide != REF_BUS) { - farr[locf + 1] += - -Bff * Vmf * Vmf + - Vmf * Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); - flps += 9 + EXAGO_FLOPS_SINOP + EXAGO_FLOPS_COSOP; - } - } else { - if (bust->ide != REF_BUS) { - farr[loct] += Gtt * Vmt * Vmt + - Vmt * Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); - flps += 9 + EXAGO_FLOPS_SINOP + EXAGO_FLOPS_COSOP; + if (!line->isdcline) { + PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt; + Gff = line->yff[0]; + Bff = line->yff[1]; + Gft = line->yft[0]; + Bft = line->yft[1]; + Gtf = line->ytf[0]; + Btf = line->ytf[1]; + Gtt = line->ytt[0]; + Btt = line->ytt[1]; + + thetaft = thetaf - thetat; + thetatf = thetat - thetaf; + + if (bus == busf) { + if (busf->ide != REF_BUS) { + farr[locf] += Gff * Vmf * Vmf + + Vmf * Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); + flps += 9 + EXAGO_FLOPS_SINOP + EXAGO_FLOPS_COSOP; + } + if (busf->ide != PV_BUS && busf->ide != REF_BUS) { + farr[locf + 1] += + -Bff * Vmf * Vmf + + Vmf * Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); + flps += 9 + EXAGO_FLOPS_SINOP + EXAGO_FLOPS_COSOP; + } + } else { + if (bust->ide != REF_BUS) { + farr[loct] += Gtt * Vmt * Vmt + + Vmt * Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); + flps += 9 + EXAGO_FLOPS_SINOP + EXAGO_FLOPS_COSOP; + } + if (bust->ide != PV_BUS && bust->ide != REF_BUS) { + farr[loct + 1] += + -Btt * Vmt * Vmt + + Vmt * Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); + flps += 9 + EXAGO_FLOPS_SINOP + EXAGO_FLOPS_COSOP; + } } - if (bust->ide != PV_BUS && bust->ide != REF_BUS) { - farr[loct + 1] += - -Btt * Vmt * Vmt + - Vmt * Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); - flps += 9 + EXAGO_FLOPS_SINOP + EXAGO_FLOPS_COSOP; + } else if (line->isdcline) { + /* DC line */ + if (bus == busf) { /* From bus */ + farr[locf] += line->pf; + farr[locf + 1] += Vmf - line->Vf; + flps += 3; + } else { + /* To bus */ + double loss; + loss = line->loss0 + line->loss1 * line->pf; + farr[loct] -= line->pf - loss; + farr[loct + 1] += Vmt - line->Vt; + flps += 6; } } } - flps += 2 * nconnlines; // for the two flops line: 447,448 } ierr = VecRestoreArrayRead(localX, &xarr); @@ -794,7 +837,9 @@ PetscErrorCode PFLOWSolutionToPS(PFLOW pflow) { for (i = 0; i < ps->nline; i++) { line = &ps->line[i]; - line->pf = line->qf = line->pt = line->qt = 0.0; + if (!line->isdcline) { + line->pf = line->qf = line->pt = line->qt = 0.0; + } } for (i = 0; i < ps->nbus; i++) { @@ -860,15 +905,6 @@ PetscErrorCode PFLOWSolutionToPS(PFLOW pflow) { line = connlines[k]; if (!line->status) continue; - PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt; - Gff = line->yff[0]; - Bff = line->yff[1]; - Gft = line->yft[0]; - Bft = line->yft[1]; - Gtf = line->ytf[0]; - Btf = line->ytf[1]; - Gtt = line->ytt[0]; - Btt = line->ytt[1]; const PSBUS *connbuses; PSBUS busf, bust; @@ -892,22 +928,44 @@ PetscErrorCode PFLOWSolutionToPS(PFLOW pflow) { thetaft = thetaf - thetat; thetatf = thetat - thetaf; - if (bus == busf) { - line->pf = Gff * Vmf * Vmf + - Vmf * Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); - line->qf = -Bff * Vmf * Vmf + - Vmf * Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); - Pnet += line->pf; - Qnet += line->qf; - } else { - line->pt = Gtt * Vmt * Vmt + - Vmt * Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); - line->qt = -Btt * Vmt * Vmt + - Vmt * Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); - Pnet += line->pt; - Qnet += line->qt; + if (!line->isdcline) { + PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt; + Gff = line->yff[0]; + Bff = line->yff[1]; + Gft = line->yft[0]; + Bft = line->yft[1]; + Gtf = line->ytf[0]; + Btf = line->ytf[1]; + Gtt = line->ytt[0]; + Btt = line->ytt[1]; + + if (bus == busf) { + line->pf = Gff * Vmf * Vmf + + Vmf * Vmt * (Gft * cos(thetaft) + Bft * sin(thetaft)); + line->qf = -Bff * Vmf * Vmf + + Vmf * Vmt * (-Bft * cos(thetaft) + Gft * sin(thetaft)); + Pnet += line->pf; + Qnet += line->qf; + } else { + line->pt = Gtt * Vmt * Vmt + + Vmt * Vmf * (Gtf * cos(thetatf) + Btf * sin(thetatf)); + line->qt = -Btt * Vmt * Vmt + + Vmt * Vmf * (-Btf * cos(thetatf) + Gtf * sin(thetatf)); + Pnet += line->pt; + Qnet += line->qt; + } + flps += 18 + 2 * (EXAGO_FLOPS_SINOP + EXAGO_FLOPS_COSOP); + } else if (line->isdcline) { + if (bus == busf) { + Pnet += line->pf; + } else { + double loss; + loss = line->loss0 + line->loss1 * line->pf; + line->pt = (line->pf - loss); + + Pnet += line->pt; + } } - flps += 18 + 2 * (EXAGO_FLOPS_SINOP + EXAGO_FLOPS_COSOP); } flps += 2 * nconnlines; diff --git a/src/ps/ps.cpp b/src/ps/ps.cpp index 4a89bac85..1a03cba0a 100644 --- a/src/ps/ps.cpp +++ b/src/ps/ps.cpp @@ -772,8 +772,11 @@ PetscErrorCode PSCreate(MPI_Comm mpicomm, PS *psout) { ps->Ngen = -1; ps->NgenON = -1; ps->Nline = -1; + ps->Ndcline = -1; ps->nlineON = -1; + ps->ndclineON = -1; ps->NlineON = -1; + ps->NdclineON = -1; ps->Nload = -1; ps->refct = 0; ps->app = NULL; @@ -972,11 +975,25 @@ PetscErrorCode PSSetUp(PS ps) { CHKERRQ(ierr); /* Set up edge connectivity */ - + /* Edges include true edges + dummy edges set up for isolated buses + These edges for isolated buses are inserted so that DMNetwork + correctly knows the number of vertices. DMNetwork filters out + any nodes that do not have edges and so this messes up our network + data layout that has these isolated buses. So, we create dummy + edges connecting these isolated buses. These are only used here + for the DMNetwork business, not anywhere else + */ ierr = PetscCalloc1(2 * Nlines, &lineconn[0]); CHKERRQ(ierr); ierr = PSGetLineConnectivity(ps, Nlines, lineconn[0]); CHKERRQ(ierr); + /* Insert dummy edges for isolated_buses */ + // for(i=0; i < ps->nisolated_buses; i++) { + /* Each dummy line is connected between bus 0 and the isolated bus */ + // lineconn[0][2 * (Nlines+i) ] = ps->isolated_buses[i]; + // lineconn[0][2 * (Nlines+i) + 1] = 0; + // } + // ierr = PetscFree(ps->isolated_buses); /* Set sizes for the network and provide edge connectivity information */ ierr = DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, 1); @@ -1068,9 +1085,10 @@ PetscErrorCode PSSetUp(PS ps) { CHKERRQ(ierr); } - /* Broadcast global Nbus,Ngen,Nbranch, Nload,and maxbusnum to all processors + /* Broadcast global Nbus,Ngen,Nbranch, Nload,Ndcline, and maxbusnum to all + * processors */ - PetscInt temp[8]; + PetscInt temp[9]; /* Pack variables */ temp[0] = ps->Nbus; temp[1] = ps->Ngen; @@ -1080,6 +1098,7 @@ PetscErrorCode PSSetUp(PS ps) { temp[5] = ps->NgenON; temp[6] = ps->NlineON; temp[7] = ps->Nref; + temp[8] = ps->Ndcline; ierr = MPI_Bcast(temp, 8, MPI_INT, 0, ps->comm->type); CHKERRQ(ierr); /* Unpack */ @@ -1091,6 +1110,7 @@ PetscErrorCode PSSetUp(PS ps) { ps->NgenON = temp[5]; ps->NlineON = temp[6]; ps->Nref = temp[7]; + ps->Ndcline = temp[8]; /* Recreate busext2intmap..this will map the local bus numbers to external * numbers */ @@ -1099,7 +1119,8 @@ PetscErrorCode PSSetUp(PS ps) { for (i = 0; i < ps->maxbusnum + 1; i++) ps->busext2intmap[i] = -1; - ps->ngen = ps->nload = ps->ngenON = ps->nlineON = 0; + ps->ngen = ps->nload = ps->ndcline = ps->ngenON = ps->nlineON = + ps->ndclineON = 0; /* Get the local number of gens and loads */ for (i = 0; i < nv; i++) { ierr = DMNetworkGetNumComponents(ps->networkdm, vtx[i], &numComponents); @@ -1133,6 +1154,10 @@ PetscErrorCode PSSetUp(PS ps) { ierr = PetscMemcpy(&ps->line[i], component, sizeof(struct _p_PSLINE)); CHKERRQ(ierr); ps->nlineON += ps->line[i].status; + if (ps->line[i].isdcline) { + ps->ndcline++; + ps->ndclineON += ps->line[i].status; + } } PetscInt genj = 0, loadj = 0; PetscInt genctr, loadctr; @@ -1170,6 +1195,7 @@ PetscErrorCode PSSetUp(PS ps) { ps->nload = ps->Nload; ps->ngenON = ps->NgenON; ps->nlineON = ps->NlineON; + ps->ndclineON = ps->NdclineON; } /* Set up @@ -1178,8 +1204,9 @@ PetscErrorCode PSSetUp(PS ps) { (c) incident generators at bus (d) incident loads at bus (e) kv levels for lines - (e) sets the starting location for the variables for this bus in the given + (f) sets the starting location for the variables for this bus in the given application state vector + (g) marks from and to buses for DC lines (ON) as PV buses */ PetscInt eStart, eEnd, vStart, vEnd; PetscInt nlines, k; @@ -1279,6 +1306,42 @@ PetscErrorCode PSSetUp(PS ps) { /* Get KV levels */ ierr = PSGetKVLevels(ps, &ps->nkvlevels, (const PetscScalar **)&ps->kvlevels); CHKERRQ(ierr); + + PSLINE line; + /* Set up for lines + - Turn off lines that are connected to isolated buses (if they are not OFF + already) + - Mark DC line ends as PV buses */ + for (i = 0; i < ps->nline; i++) { + line = &ps->line[i]; + + const PSBUS *connbuses; + PSBUS busf, bust; + /* Get the connected buses */ + ierr = PSLINEGetConnectedBuses(line, &connbuses); + CHKERRQ(ierr); + busf = connbuses[0]; + bust = connbuses[1]; + + if (busf->ide == ISOLATED_BUS || bust->ide == ISOLATED_BUS) { + /* Switch off lines connected to isolated buses */ + if (line->status) { + line->status = 0; + ps->nlineON--; + } + } + + if (!line->isdcline) + continue; + if (!line->status) + continue; + + if (busf->ide != REF_BUS || busf->ide != PV_BUS) + busf->ide = PV_BUS; + if (bust->ide != REF_BUS || bust->ide != PV_BUS) + bust->ide = PV_BUS; + } + // ierr = PetscPrintf(PETSC_COMM_SELF,"Rank %d Came // here\n",ps->comm->rank);CHKERRQ(ierr); ps->setupcalled = PETSC_TRUE; diff --git a/src/ps/psreaddata.cpp b/src/ps/psreaddata.cpp index 5c4405b5d..2f6f7bf1a 100644 --- a/src/ps/psreaddata.cpp +++ b/src/ps/psreaddata.cpp @@ -136,7 +136,7 @@ PetscErrorCode PSReadPSSERawData(PS ps, const char netfile[]) { ps->Ngen = ps->ngen = Ngenerator; ps->Nline = ps->nline = Nline + Ntransformer; ps->NgenON = 0; - ps->NlineON = 0; + ps->nlineON = ps->NlineON = 0; #if defined DEBUGPS ierr = PetscPrintf( @@ -480,6 +480,20 @@ PetscErrorCode PSReadPSSERawData(PS ps, const char netfile[]) { PetscFunctionReturn(0); } +/* Returns true if busum is in the isolated_buses list */ +PetscBool BusIsIsolated(PetscInt busnum, PetscInt nisolated_buses, + PetscInt *isolated_buses) { + PetscInt i; + + for (i = 0; i < nisolated_buses; i++) { + if (isolated_buses[i] == busnum) { + return PETSC_TRUE; + } + } + + return PETSC_FALSE; +} + /* PSReadMatPowerData - Reads the MATPOWER data file and populates the PS object @@ -502,13 +516,20 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { -1; /* xx_end_line points to the next line after the record ends */ PetscInt gen_start_line = -1, gen_end_line = -1; PetscInt br_start_line = -1, br_end_line = -1; + PetscInt dcline_start_line = -1, dcline_end_line = -1; PetscInt gencost_start_line = -1, gencost_end_line = -1; PetscInt genfuel_start_line = -1, genfuel_end_line = -1; PetscInt loadcost_start_line = -1, loadcost_end_line = -1; + + /* Number of lines for each data */ + PetscInt bus_numlines, gen_numlines, br_numlines, dcline_numlines; + PetscInt load_numlines = 0; /* Number of blank lines in bus, gen, br, gencost, and genfuel branch arrays */ PetscInt bus_nblank_lines = 0, gen_nblank_lines = 0; PetscInt br_nblank_lines = 0; + PetscInt dcline_nblank_lines = 0; + char line[MAXLINE]; PetscInt loadi = 0, geni = 0, bri = 0, busi = 0, gencosti = 0, genfueli = 0, loadcosti = 0, i; @@ -517,14 +538,25 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { PetscInt intbusnum; char *str; char *out; + PetscInt buff_num = 1; + PetscInt + *temp_buff; /* buffer to hold 10 isolated buses, gets resized later */ + PetscBool bad_data = PETSC_FALSE; + const char bad_data_str[] = "Input Data Error: "; + PetscInt num_max_cost_coeffs = + 3; // Max. number of cost-coefficients for any generator PetscFunctionBegin; if (ps->comm->type != PETSC_COMM_SELF && ps->comm->rank != 0) { - ps->Nline = ps->Nbus = ps->Ngen = ps->Nload = 0; + ps->Nline = ps->Nbus = ps->Ngen = ps->Nload = ps->Ndcline = 0; PetscFunctionReturn(0); } + ps->nisolated_buses = 0; + ierr = PetscCalloc1(10, &ps->isolated_buses); + CHKERRQ(ierr); + fp = fopen(netfile, "r"); /* Check for valid file */ if (fp == NULL) { @@ -554,6 +586,9 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { line_counter + 1; /* Generator data starts from next line */ if (strstr(line, "mpc.branch") != NULL) br_start_line = line_counter + 1; /* Branch data starts from next line */ + if (strstr(line, "mpc.dcline") != NULL && dcline_start_line == -1) + dcline_start_line = + line_counter + 1; /* DC line data starts from next line */ if (strstr(line, "mpc.gencost") != NULL) gencost_start_line = line_counter + 1; /* Gen cost data starts from next line */ @@ -580,6 +615,8 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { gen_end_line = line_counter; if (br_start_line != -1 && br_end_line == -1) br_end_line = line_counter; + if (dcline_start_line != -1 && dcline_end_line == -1) + dcline_end_line = line_counter; if (gencost_start_line != -1 && gencost_end_line == -1) gencost_end_line = line_counter; if (loadcost_start_line != -1 && loadcost_end_line == -1) @@ -588,59 +625,90 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { if (bus_start_line != -1 && bus_end_line == -1) { if (strcmp(line, "\n") == 0 || strcmp(line, "\r\n") == 0) - bus_nblank_lines++; + bus_nblank_lines = bus_nblank_lines + 1; } if (gen_start_line != -1 && gen_end_line == -1) { if (strcmp(line, "\n") == 0 || strcmp(line, "\r\n") == 0) - gen_nblank_lines++; + gen_nblank_lines = gen_nblank_lines + 1; } if (br_start_line != -1 && br_end_line == -1) { if (strcmp(line, "\n") == 0 || strcmp(line, "\r\n") == 0) - br_nblank_lines++; + br_nblank_lines = br_nblank_lines + 1; + } + + if (dcline_start_line != -1 && dcline_end_line == -1) { + if (strcmp(line, "\n") == 0 || strcmp(line, "\r\n") == 0) + dcline_nblank_lines = dcline_nblank_lines + 1; } /* Count the number of pq loads */ if (bus_start_line != -1 && line_counter >= bus_start_line && bus_end_line == -1) { sscanf(line, "%d %d %lf %lf", &extbusnum, &bustype_i, &Pd, &Qd); + if (bustype_i == ISOLATED_BUS) { /* Isolated bus */ + /* put it in the isolated_buses array for later work */ + ps->isolated_buses[ps->nisolated_buses++] = extbusnum; + if (ps->nisolated_buses == buff_num * 10) { + /* Reached buffer max, need to allocate more size */ + buff_num++; + ierr = PetscCalloc1(buff_num * 10, &temp_buff); + CHKERRQ(ierr); + /* Copy over isolated_buses to temp_buff */ + ierr = PetscMemcpy(temp_buff, ps->isolated_buses, + (buff_num - 1) * 10 * sizeof(PetscInt)); + ierr = PetscFree(ps->isolated_buses); + CHKERRQ(ierr); + /* Updated pointer to isolated_buses */ + ps->isolated_buses = temp_buff; + } + } + if (!((Pd == 0.0) && (Qd == 0.0))) - ps->Nload++; + load_numlines++; if (extbusnum > ps->maxbusnum) ps->maxbusnum = extbusnum; } + + /* Count the number of cost-coefficients */ + if (gencost_start_line != -1 && line_counter >= gencost_start_line && + gencost_end_line == -1) { + PetscInt gencost_model, num_cost_coeffs; + PetscScalar cost_startup, cost_shutdown; + sscanf(line, "%d %lf %lf %d", &gencost_model, &cost_startup, + &cost_shutdown, &num_cost_coeffs); + if (num_cost_coeffs > num_max_cost_coeffs) + num_max_cost_coeffs = num_cost_coeffs; + } + line_counter++; } fclose(fp); - ps->Nbus = ps->nbus = bus_end_line - bus_start_line - bus_nblank_lines; - ps->Ngen = ps->ngen = gen_end_line - gen_start_line - gen_nblank_lines; - ps->Nline = ps->nline = br_end_line - br_start_line - br_nblank_lines; - ps->nload = ps->Nload; + bus_numlines = bus_end_line - bus_start_line - bus_nblank_lines; + gen_numlines = gen_end_line - gen_start_line - gen_nblank_lines; + br_numlines = br_end_line - br_start_line - br_nblank_lines; + dcline_numlines = dcline_end_line - dcline_start_line - dcline_nblank_lines; + ps->NgenON = 0; - ps->NlineON = 0; + ps->nlineON = ps->NlineON = 0; + ps->ndclineON = ps->NdclineON = 0; -#if defined DEBUGPS - ierr = PetscPrintf( - PETSC_COMM_SELF, - "System summary : Nbuses = %d, Ngen = %d, Nload = %d, Nbranch = %d\n", - ps->Nbus, ps->Ngen, ps->Nload, ps->Nline); + ierr = PetscCalloc1(bus_numlines, &ps->bus); CHKERRQ(ierr); -#endif - ierr = PetscCalloc1(ps->Nbus, &ps->bus); + ierr = PetscCalloc1(gen_numlines, &ps->gen); CHKERRQ(ierr); - ierr = PetscCalloc1(ps->Ngen, &ps->gen); + ierr = PetscCalloc1(load_numlines, &ps->load); CHKERRQ(ierr); - ierr = PetscCalloc1(ps->Nload, &ps->load); - CHKERRQ(ierr); - ierr = PetscCalloc1(ps->Nline, &ps->line); + ierr = PetscCalloc1(br_numlines + dcline_numlines, + &ps->line); /* Includes space for DC lines */ CHKERRQ(ierr); Bus = ps->bus; Gen = ps->gen; Load = ps->load; Branch = ps->line; - for (i = 0; i < ps->Nbus; i++) { + for (i = 0; i < bus_numlines; i++) { ps->bus[i].ngen = ps->bus[i].nload = ps->bus[i].ngenON = ps->bus[i].nshunt = 0; ps->bus[i].qrange = ps->bus[i].qmintot = ps->bus[i].Pgtot = @@ -658,6 +726,23 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { for (i = 0; i < ps->maxbusnum + 1; i++) busext2intmap[i] = -1; + ps->Nbus = ps->nbus = 0; + ps->Nline = ps->nline = 0; + ps->Nload = ps->nload = 0; + ps->Ngen = ps->ngen = 0; + ps->Ndcline = ps->ndcline = 0; + + /* Create a temp array to store the locations in mpc.gen (generators) + connected to isolated buses. This is used later in ignoring these + generators in gencost data + */ + PetscInt + temp_geni[200]; // Assume there are max 200 generators at an isolated node + for (i = 0; i < 200; i++) + temp_geni[i] = line_counter + 1; // initialize + + PetscInt i_idx = -1, isol_idx = 0, isol_idx2 = 0, isol_idx3 = 0; + fp = fopen(netfile, "r"); /* Reading data */ for (i = 0; i < line_counter; i++) { @@ -673,8 +758,23 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { &Bus[busi].basekV, &Bus[busi].zone, &Bus[busi].Vmax, &Bus[busi].Vmin); + if (Bus[busi].ide == ISOLATED_BUS) + continue; /* Skip isolated bus */ + + ps->nbus++; + ps->Nbus++; + Bus[busi].Vmax = Bus[busi].Vmax == 0 ? 1.1 : Bus[busi].Vmax; Bus[busi].Vmin = Bus[busi].Vmin == 0 ? 0.9 : Bus[busi].Vmin; + /* Sanity check for voltage limits */ + if ((Bus[busi].Vmax < Bus[busi].Vmin) || (Bus[busi].Vmax > 1.1) || + (Bus[busi].Vmin < 0.9)) { + bad_data = PETSC_TRUE; + ierr = PetscPrintf( + ps->comm->type, + "%sVmax = %4.3f at bus %d is less than Vmin = %4.3f\n", + bad_data_str, Bus[busi].Vmax, Bus[busi].bus_i, Bus[busi].Vmin); + } // Convert angle to radians Bus[busi].va *= PETSC_PI / 180.0; @@ -689,6 +789,9 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { Bus[busi].nshunt++; if (!((Pd == 0.0) && (Qd == 0.0))) { + ps->nload++; + ps->Nload++; + Load[loadi].bus_i = Bus[busi].bus_i; Load[loadi].status = 1; Load[loadi].pl = Pd / ps->MVAbase; @@ -710,13 +813,16 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { "Exceeded maximum number of loads allowed at bus"); loadi++; } + Bus[busi].nconnlines = 0; busi++; } /* Read generator data */ if (i >= gen_start_line && i < gen_end_line) { + i_idx += 1; if (strcmp(line, "\n") == 0 || strcmp(line, "\r\n") == 0) continue; + sscanf(line, "%d %lf %lf %lf %lf %lf %lf %d %lf %lf %lf %lf %lf %lf %lf %lf " "%lf %lf %lf %lf %lf", @@ -728,6 +834,19 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { &Gen[geni].ramp_rate_10min, &Gen[geni].ramp_rate_30min, &Gen[geni].ramp_rate_min_mvar, &Gen[geni].apf); + if (BusIsIsolated(Gen[geni].bus_i, ps->nisolated_buses, + ps->isolated_buses)) { + if (isol_idx == 200) { + SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, + "Total number of generators at isolated buses greater than " + "max. allowed = 200\n"); + } + temp_geni[isol_idx++] = i_idx; + continue; /* Skip generator at isolated bus */ + } + ps->ngen++; + ps->Ngen++; + intbusnum = busext2intmap[Gen[geni].bus_i]; Gen[geni].qt = Gen[geni].qt > 1e10 ? PETSC_INFINITY : Gen[geni].qt; @@ -767,8 +886,8 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { Gen[geni].genfuel_type = GENFUEL_UNDEFINED; Gen[geni].ramp_rate_min = GENRAMPRATE_COAL / ps->MVAbase; /* Defaults to COAL ramp rate */ - Gen[geni].ramp_rate_10min = Gen[genfueli].ramp_rate_min * 10; - Gen[geni].ramp_rate_30min = Gen[genfueli].ramp_rate_min * 30; + Gen[geni].ramp_rate_10min = Gen[geni].ramp_rate_min * 10; + Gen[geni].ramp_rate_30min = Gen[geni].ramp_rate_min * 30; } if (Gen[geni].status) { @@ -788,6 +907,26 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { */ snprintf(Gen[geni].id, 3, "%-2d", 1 + Bus[intbusnum].ngen); + /* Sanity check for real power limits */ + if (Gen[geni].pt < Gen[geni].pb) { + bad_data = PETSC_TRUE; + ierr = PetscPrintf(ps->comm->type, + "%sPmax = %4.3f for generator with id %s at bus %d " + "is less than Pmin = %4.3f\n", + bad_data_str, Gen[geni].pt, Gen[geni].id, + Gen[geni].bus_i, Gen[geni].pb); + } + + /* Sanity check for reactive power limits */ + if (Gen[geni].qt < Gen[geni].qb) { + bad_data = PETSC_TRUE; + ierr = PetscPrintf(ps->comm->type, + "%sQmax = %4.3f for generator with id %s at bus %d " + "is less than Qmin = %4.3f\n", + bad_data_str, Gen[geni].qt, Gen[geni].id, + Gen[geni].bus_i, Gen[geni].qb); + } + Bus[intbusnum].gidx[Bus[intbusnum].ngen++] = geni; // Bus[intbusnum].vm = Gen[geni].vs; @@ -809,10 +948,51 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { if (i >= gencost_start_line && i < gencost_end_line) { if (strcmp(line, "\n") == 0 || strcmp(line, "\r\n") == 0) continue; - sscanf(line, "%d %lf %lf %d %lf %lf %lf", &Gen[gencosti].cost_model, - &Gen[gencosti].cost_startup, &Gen[gencosti].cost_shutdown, - &Gen[gencosti].cost_ncoeffs, &Gen[gencosti].cost_alpha, - &Gen[gencosti].cost_beta, &Gen[gencosti].cost_gamma); + + if (temp_geni[isol_idx2] == i - gencost_start_line) { + isol_idx2 += 1; + continue; + } + + if (num_max_cost_coeffs == 2) { + Gen[geni].cost_alpha = 0.0; + sscanf(line, "%d %lf %lf %d %lf %lf", &Gen[gencosti].cost_model, + &Gen[gencosti].cost_startup, &Gen[gencosti].cost_shutdown, + &Gen[gencosti].cost_ncoeffs, &Gen[gencosti].cost_beta, + &Gen[gencosti].cost_gamma); + } else if (num_max_cost_coeffs == 3) { + sscanf(line, "%d %lf %lf %d %lf %lf %lf", &Gen[gencosti].cost_model, + &Gen[gencosti].cost_startup, &Gen[gencosti].cost_shutdown, + &Gen[gencosti].cost_ncoeffs, &Gen[gencosti].cost_alpha, + &Gen[gencosti].cost_beta, &Gen[gencosti].cost_gamma); + } else if (num_max_cost_coeffs == 4) { + sscanf(line, "%d %lf %lf %d %*f %lf %lf %lf", &Gen[gencosti].cost_model, + &Gen[gencosti].cost_startup, &Gen[gencosti].cost_shutdown, + &Gen[gencosti].cost_ncoeffs, &Gen[gencosti].cost_alpha, + &Gen[gencosti].cost_beta, &Gen[gencosti].cost_gamma); + } else { + bad_data = PETSC_TRUE; + ierr = + PetscPrintf(ps->comm->type, + "%sExaGO can only read generator cost curves upto 3rd " + "order. More than 4 coefficients found for cost curve " + "of cost curve for generator with id %s at bus %d\n", + bad_data_str, Gen[gencosti].id, Gen[gencosti].bus_i); + CHKERRQ(ierr); + } + + /* Sanity check for generator cost function */ + if (Gen[gencosti].cost_model != 2) { + bad_data = PETSC_TRUE; + ierr = + PetscPrintf(ps->comm->type, + "%sExaGO only supports generator cost curves modeled " + "as quadratic function.. No quadratic cost curve " + "function found for generator with id %s at bus %d\n", + bad_data_str, Gen[gencosti].id, Gen[gencosti].bus_i); + CHKERRQ(ierr); + } + gencosti++; } @@ -830,6 +1010,10 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { if (i >= genfuel_start_line && i < genfuel_end_line) { if (strcmp(line, "\n") == 0 || strcmp(line, "\r\n") == 0) continue; + if (temp_geni[isol_idx3] == i - genfuel_start_line) { + isol_idx3 += 1; + continue; + } if (strstr(line, "coal") != NULL) { Gen[genfueli].genfuel_type = GENFUEL_COAL; Gen[genfueli].ramp_rate_min = GENRAMPRATE_COAL / ps->MVAbase; @@ -841,8 +1025,6 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { Gen[genfueli].ramp_rate_min = GENRAMPRATE_WIND / ps->MVAbase; Gen[genfueli].ramp_rate_10min = Gen[genfueli].ramp_rate_min * 10; Gen[genfueli].ramp_rate_30min = Gen[genfueli].ramp_rate_min * 30; - Gen[genfueli].pb = 0.0; /* Set lower Pg limit to 0.0 so that wind power - can be curtailed if need be */ ps->ngenwind++; ps->ngenrenew++; Gen[genfueli].isrenewable = PETSC_TRUE; @@ -895,8 +1077,20 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { &Branch[bri].x, &Branch[bri].b, &Branch[bri].rateA, &Branch[bri].rateB, &Branch[bri].rateC, &Branch[bri].tapratio, &Branch[bri].phaseshift, &Branch[bri].status); - if (Branch[bri].status) + + if ((BusIsIsolated(Branch[bri].fbus, ps->nisolated_buses, + ps->isolated_buses)) || + (BusIsIsolated(Branch[bri].tbus, ps->nisolated_buses, + ps->isolated_buses))) + continue; /* Skip branch connected to isolated bus */ + + ps->nline++; + ps->Nline++; + + if (Branch[bri].status) { + ps->nlineON++; ps->NlineON++; + } if (!Branch[bri].tapratio) Branch[bri].tapratio = 1.0; Branch[bri].phaseshift *= PETSC_PI / 180.0; @@ -922,6 +1116,10 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { intbusnum = busext2intmap[Branch[bri].tbus]; Branch[bri].internal_j = intbusnum; + /* Update number of connected lines for the from and to bus */ + Bus[Branch[bri].internal_i].nconnlines++; + Bus[Branch[bri].internal_j].nconnlines++; + PetscInt lineididx = 0; for (linenum = 0; linenum < bri - 1; linenum++) { if (Branch[bri].internal_i == Branch[linenum].internal_i && @@ -932,6 +1130,16 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { /* MatPower does not have ids for lines. Using bri+1 as the id */ snprintf(Branch[bri].ckt, 3, "%-2d", 1 + lineididx); + if (Branch[bri].rateA < 0.0) { + bad_data = PETSC_TRUE; + ierr = PetscPrintf( + ps->comm->type, + "%srateA = %4.3f for line %d -- %d with id %s is negative\n", + bad_data_str, Branch[bri].rateA, Branch[bri].fbus, Branch[bri].tbus, + Branch[bri].ckt); + CHKERRQ(ierr); + } + /* Compute self and transfer admittances */ PetscScalar R, X, Bc, B, G, Zm, tap, shift, tap2, tapr, tapi; R = Branch[bri].r; @@ -981,9 +1189,111 @@ PetscErrorCode PSReadMatPowerData(PS ps, const char netfile[]) { Branch[bri].bdc = (1 / X) / tap; Branch[bri].pshift = -(1 / X) * shift; } + Branch[bri].isdcline = PETSC_FALSE; + bri++; + } + + /* DC lines */ + if (i >= dcline_start_line && i < dcline_end_line) { + if (strcmp(line, "\n") == 0 || strcmp(line, "\r\n") == 0) + continue; + sscanf(line, + "%d %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", + &Branch[bri].fbus, &Branch[bri].tbus, &Branch[bri].status, + &Branch[bri].pf, &Branch[bri].pt, &Branch[bri].qf, &Branch[bri].qt, + &Branch[bri].Vf, &Branch[bri].Vt, &Branch[bri].pmin, + &Branch[bri].pmax, &Branch[bri].qminf, &Branch[bri].qmaxf, + &Branch[bri].qmint, &Branch[bri].qmaxt, &Branch[bri].loss0, + &Branch[bri].loss1); + + if ((BusIsIsolated(Branch[bri].fbus, ps->nisolated_buses, + ps->isolated_buses)) || + (BusIsIsolated(Branch[bri].tbus, ps->nisolated_buses, + ps->isolated_buses))) + continue; /* Skip branch connected to isolated bus */ + ps->ndcline++; + ps->Ndcline++; + + Branch[bri].pf /= ps->MVAbase; + Branch[bri].pmin /= ps->MVAbase; + Branch[bri].pmax /= ps->MVAbase; + Branch[bri].qminf /= ps->MVAbase; + Branch[bri].qmaxf /= ps->MVAbase; + Branch[bri].qmint /= ps->MVAbase; + Branch[bri].qmaxt /= ps->MVAbase; + Branch[bri].loss0 /= ps->MVAbase; + + Branch[bri].isdcline = PETSC_TRUE; + ps->Nline++; + ps->nline++; + if (Branch[bri].status) { + ps->NlineON++; + ps->nlineON++; + ps->NdclineON++; + ps->ndclineON++; + } + + intbusnum = busext2intmap[Branch[bri].fbus]; + Branch[bri].internal_i = intbusnum; + + intbusnum = busext2intmap[Branch[bri].tbus]; + Branch[bri].internal_j = intbusnum; + + /* Update number of connected lines for the from and to bus */ + Bus[Branch[bri].internal_i].nconnlines++; + Bus[Branch[bri].internal_j].nconnlines++; + + PetscInt lineididx = 0; + for (linenum = 0; linenum < bri - 1; linenum++) { + if (Branch[bri].internal_i == Branch[linenum].internal_i && + Branch[bri].internal_j == Branch[linenum].internal_j) + lineididx += 1; + } + + /* MatPower does not have ids for lines. Using bri+1 as the id */ + snprintf(Branch[bri].ckt, 3, "%-2d", 1 + lineididx); + + if (Branch[bri].pmax < Branch[bri].pmin) { + bad_data = PETSC_TRUE; + ierr = PetscPrintf(ps->comm->type, + "%sPmax = %4.3f for DC line %d -- %d with id %s is " + "less than Pmin = %4.3f\n", + bad_data_str, Branch[bri].pmax, Branch[bri].fbus, + Branch[bri].tbus, Branch[bri].ckt, Branch[bri].pmin); + CHKERRQ(ierr); + } + + if (Branch[bri].qmaxf < Branch[bri].qminf) { + bad_data = PETSC_TRUE; + ierr = + PetscPrintf(ps->comm->type, + "%sQmaxf = %4.3f for DC line %d -- %d with id %s is " + "less than Qminf = %4.3f\n", + bad_data_str, Branch[bri].qmaxf, Branch[bri].fbus, + Branch[bri].tbus, Branch[bri].ckt, Branch[bri].qminf); + CHKERRQ(ierr); + } + + if (Branch[bri].qmaxt < Branch[bri].qmint) { + bad_data = PETSC_TRUE; + ierr = + PetscPrintf(ps->comm->type, + "%sQmaxt = %4.3f for DC line %d -- %d with id %s is " + "less than Qmint = %4.3f\n", + bad_data_str, Branch[bri].qmaxt, Branch[bri].fbus, + Branch[bri].tbus, Branch[bri].ckt, Branch[bri].qmint); + CHKERRQ(ierr); + } + bri++; } } + + if (bad_data) { + SETERRQ(ps->comm->type, 0, + "Found errors in input data. Please fix the suggested errors"); + } + fclose(fp); PetscFunctionReturn(0); diff --git a/tests/functionality/opflow/CMakeLists.txt b/tests/functionality/opflow/CMakeLists.txt index 70a3c445c..fc54fae77 100644 --- a/tests/functionality/opflow/CMakeLists.txt +++ b/tests/functionality/opflow/CMakeLists.txt @@ -163,6 +163,22 @@ if(EXAGO_INSTALL_TESTS) HIOP ) + exago_add_test( + NAME + OPFLOW_DCLINE_TEST + COMMAND + ${RUNCMD} + $ + -netfile + ${CMAKE_SOURCE_DIR}/datafiles/case9/case9_dcline.m + -opflow_solver + IPOPT + -opflow_model + POWER_BALANCE_POLAR + DEPENDS + IPOPT + ) + foreach(fmt "MATPOWER" "CSV" "JSON" "MINIMAL") set(tname "OPFLOW_SOLUTION_OUTPUT_${fmt}") exago_add_test( diff --git a/tests/functionality/opflow/ipopt_pbpol.toml b/tests/functionality/opflow/ipopt_pbpol.toml index 375fc9a6a..2ff3ac547 100644 --- a/tests/functionality/opflow/ipopt_pbpol.toml +++ b/tests/functionality/opflow/ipopt_pbpol.toml @@ -272,3 +272,11 @@ is_powerimb_active = true gen_bus_voltage_type = 'VARIABLE_WITHIN_BOUNDS' num_iters = 27 obj_value = 740032.301353 + +[[testcase]] +network = 'datafiles/case9/case9_dcline.m' +description = "Test for OPFLOW with DC line" +is_loadloss_active = false +is_powerimb_active = false +num_iters = 18 +obj_value = 5312.16147