AC optimal power flow - #162
Open
erikfilias wants to merge 64 commits into
Open
Conversation
Adds a branch flow AC network formulation, switched on with the IndACPowerFlow option. The DC path is unchanged: every new constraint, variable and writer is behind that flag, and eBalanceElec keeps its name and index so the ten places that read its dual by string still work. What it models - Squared voltage magnitude and angle at each bus, squared current on each branch, and active and reactive power at both ends of every AC branch. - Transformer tap ratios, line charging, bus shunt devices, generator reactive capability, and synchronous condensers. - Reactive power not served, priced like energy not served, so a system short of reactive power says where instead of failing to solve. The current definition can be supplied three ways, chosen with IndACModelType: a second-order cone (0), a piecewise linear staircase (1), or the exact non-linear equation (2). The cone is a relaxation, so every run reports how far each branch sits from the boundary and warns when the gap is large. Those numbers decide whether the reported currents, losses and voltages mean anything. Bound tightening runs before the variables are declared. It narrows the voltage and angle bounds using only inequalities the model already implies, because a tighter number that the data does not support would cut off the true optimum. New files - openTEPES_InputDataAC.py, openTEPES_BoundTightening.py, openTEPES_SettingUpVariablesAC.py, openTEPES_ModelFormulationAC.py and openTEPES_OutputResultsAC.py. - Two cases, 9n_AC and RTS-GMLC_AC. - tests/test_ac_input.py, 36 tests. - Design notes under doc/design and the prototypes the choices were measured on. Combinations that are refused rather than solved wrongly: cycle flow, single node, variable TTC and PTDF all conflict with the AC model, and each raises a clear error. Known limitations - RTS-GMLC_AC has not been verified end to end. The case is a full year at hourly resolution, 8736 load levels over 73 buses and 120 branches, which builds a model of 27.5 million rows and 40 million columns. A barrier solve was still running after nine hours, which is unremarkable at that size. The matrix range is 1e-05 to 3e+02, so the model is well conditioned and the cost is size rather than numerics. The AC formulation writes 18 constraints per branch per hour where the DC one writes 7, on top of six more variables per branch. - The AC current penalty is 1e-3 and is not a tie-breaker at that size. On 9n_AC it is 2.4 per cent of total system cost and it enters the prices reported as locational marginal costs. It has its own row in the cost summary so it can be seen. Section 11 of doc/design/AC_OPF_Prototype_Results.md has the measurements and the trade-off. - Validation against an exact AC power flow has not been run. It needs ipopt. - Shunt and condenser candidates are not supported by the Benders decomposition path, which raises an error rather than building a wrong master problem.
RTS-GMLC_AC carries no candidates of any kind, so it is already an operation
only case. What makes it a 27.5 million row model is the 8736 hour horizon.
A barrier solve had not finished after nine hours, which makes it unusable for
day to day work on the AC model.
These two cases are the same 168 hours over the same network, one under each
network model, so the cost of AC can be read directly:
rows columns nonzeros wall
DC 98,343 123,227 423,978 8.9 s
AC 528,423 770,531 2,424,390 103.7 s
AC / DC 5.4x 6.3x 5.7x 11.7x
The week is the one containing the annual peak, 23 to 30 August, peaking at
10,212 MW. Time indexed tables are trimmed rather than zeroed, so each case is
about 0.9 MB against 45 MB for the full year. Weights are left at 1: the costs
these cases report are for the week they model, not an annual figure.
The AC case also shows something the nine bus case cannot. The conic relaxation
is not tight on 21 of its 120 branches, and on three of them the slack is most
of the branch rating. The current, loss and loading reported on those branches
are not supported by the flows and should not be read as physical. The three
worst are the same branch position in each of the three RTS areas, which points
at the network structure rather than at a data error. Section 12 of
doc/design/AC_OPF_Prototype_Results.md has the numbers.
This is the case for keeping the relaxation diagnostic in the minimal output
set. Without it a user reading only the voltages and currents has no way to know
that a sixth of the branches carry numbers that do not mean what they look like.
Two changes. The first is a bug in what was already committed. The angle relation had the wrong sign eAngleEnvM built the envelope numerator as x*P + r*Q. It is a minus. From S_ij = V_i' conj((V_i' - V_j) y) with y = (r - jx)/z^2: P = [(v_i^2 - v_i v_j cos th) r + (v_i v_j sin th) x] / z^2 Q = [(v_i^2 - v_i v_j cos th) x - (v_i v_j sin th) r] / z^2 so x P - r Q = v_i v_j sin th. Measured against the textbook pi-model computed from the model's own voltages and angles, the worst branch flow error goes from 38 MW to 0.00001 MW, and the angles now close around every cycle to 6e-18 radians instead of 0.634 degrees. Before the fix no assignment of bus angles reproduced the flows, so the reported operating point was not one the system could take up. Nothing caught this for ten review rounds because the relaxation gap, the angle band, the envelope and the restoration all measure the model against its own definition of that relation. The wrong sign was repeated in two module docstrings, so the comments agreed with the code. A self-consistency check cannot see a consistently wrong premise. The exact restoration pass ACRestorationPass, switched on with IndACRestore = 1. After the ordinary solve it holds the plan - commitment, switching and every investment stay where the relaxed solve put them - swaps the relaxed current definition for the exact equality and the angle envelope for the exact angle relation, and re-solves the network on ipopt. Discrete variables are fixed by domain rather than by name, so the pass stays safe when a variable is added or renamed, and any name it expects and cannot find is reported. The duals are dropped afterwards. They belong to the relaxed solve and would otherwise be published beside primal values from a different operating point. Measured on a 24 hour RTS-GMLC window, where the cone is loose on 9 of 120 branches: the relaxed solve reports 15.37 MEUR for a schedule that actually costs 24.11 once its physics are exact. It understates its own plan by 57 per cent. On 9n_AC, where the cone is tight, the restoration changes nothing and the relaxed and exact optima agree to seven significant figures. Also from the tenth review: discrete variables are held by domain, bound tightening no longer propagates across branches outside their commissioning window, a generator derated to nothing no longer supplies free reactive power, and the angle release now clears a one-sided band. Known limitations - The exact model no longer solves on RTS from a cold start: ipopt runs for 730 seconds and stops at infeasible. There is no measured exact optimum for that case. The restoration converges on the same case from the relaxed solution, which suggests a starting point rather than a genuinely empty feasible set, but the two readings have not been separated. - A Newton-Raphson comparison in pandapower still puts bus voltages 0.0130 per unit apart. The pi-model check agrees to 1e-5 MW, so the branch equations are right and the gap is in how that validation script builds the network. It is unresolved and the figure should not be quoted either way. - HVDC links are still modelled as controllable active flow with a linear loss factor. There is no converter: no reactive draw for a line-commutated station, no capability curve for a voltage-source one.
An HVDC link has a converter station at each end, and openTEPES has always
modelled the link as active power transfer with a linear loss factor and nothing
else. That is only harmless while no case has a built DC link. The two
technologies get the reactive side wrong in opposite directions, so the error
cannot be folded into the loss factor:
IndACConverter = 0 no converter, the behaviour up to now
= 1 line-commutated. Each station DRAWS reactive power at both
terminals, tan(acos(pf)) times the active power it carries,
so the AC system needs more compensation, not less.
= 2 voltage-source. Each station is a controllable reactive
source or sink within its rating, so it behaves like a
STATCOM and relieves the AC system instead.
ConverterPF sets the power factor, 0.85 by default, which gives 0.62 per unit of
transferred power. That is the usual range for a classic HVDC station.
The LCC draw needs |P_dc|, which needs the flow split in two. Pinning only the
difference is not enough, and the first version of this claimed the resulting
error was one-signed and therefore safe. It is not. The draw enters the reactive
balance with a MINUS, so at a node whose line charging exceeds its demand the
model gains by inflating both halves: it absorbs the surplus for free rather than
paying for the reactive slack, the converter becomes an unbounded reactive sink,
and a real reactive over-supply disappears from the results instead of being
reported. A direction binary per link per load level makes the split exact. Only
the LCC model pays for it, and DC links are few.
Also from the eleventh review:
- A VSC terminal is now gated on the line being in service. Without that an
unbuilt HVDC candidate handed the system a free STATCOM at both ends, so the
model would decline to build shunts and condensers it genuinely needed.
- The restored angle relation releases an out-of-service branch, as the envelope
it replaces always did. A bare equality on an unbuilt candidate forced the two
buses it would have joined to the same angle.
- The reactive marginal is no longer dropped at a node whose only connection is
an HVDC link.
- The restoration no longer warns about variables that are absent by design. It
was telling the user on nearly every run that the plan was free to move, which
was false.
- The module docstring still stated the angle relation with the sign that was
fixed in the previous commit.
Neither bundled case exercises any of this by default: 9n_AC ships a DC candidate
that is not built, and RTS-GMLC has no HVDC at all. The tests force a link into
service, and one test forbids building it and checks the converter then supplies
nothing.
…sewhere
Every other check in this repository compares the model against itself. The
relaxation gap says the cone is closed; the envelope, the band and the
restoration all measure the model against its own definition of the relations
they enforce. None of that can detect a consistently wrong premise, and one had
gone undetected through ten code reviews: the angle relation was written
x*P + r*Q when it is x*P - r*Q. It was found with this script.
prototypes/ac_formulations/validate.py runs three checks on any solved AC case,
cheapest first:
branch residual the model's own branch flows against the textbook pi-model
computed from its own voltages and angles, derived from
scratch so it depends on no openTEPES constraint
loop residual each branch's angle recovered from its own flows, summed
around every independent cycle. A solution whose loops do not
close is not an operating point the system can take up,
however tight the cone
power flow error the network rebuilt in pandapower from the same r, x, b and
tap data, the setpoints openTEPES chose injected, and
Newton-Raphson run
Results, both bundled AC cases, restoration on:
9n_AC dP 0.00001 MW loop 4e-10 rad dV 2e-09 p.u. dAngle 0
RTS-GMLC dP 0.00000 MW loop 8e-14 rad dV 2e-09 p.u. dAngle 1e-09
The second matters for the transformer tap. 9n_AC has none, so until now the tap
convention rested on a derivation and on a unit test that asserts the same
convention it was written from, which could not have caught an inverted tap.
RTS-GMLC has sixteen at 1.015 and 1.03, and pandapower's own transformer model
agrees.
The Newton-Raphson gap that the docs recorded as open was in the validation
script, not the model. It injected a candidate shunt's nameplate rating for a
device the model had declined to build, putting 299 Mvar into the network that
the solution never contained. create_impedance was suspected and was not the
cause: rebuilding the lines from explicit ohms changed the answer by nothing.
Two traps are written into the script's docstring so they are not rediscovered:
use the susceptance actually in service rather than the nameplate, and remember
that HVDC links carry active power that has to appear at both ends. A third is
recorded in the docs: zeroing the line charging halved the voltage gap and
pointed hard at the charging model, when it was the shunt error interacting with
it. An intermediate measurement that moves in the expected direction is not
evidence of the expected cause.
Not covered: both cases run here are single-period.
…laim Prose pass on this branch's own additions, per the house style: module docstrings cut from 49 and 34 lines to 17, the longest comment blocks condensed, and the three changelog entries this work was missing. Also corrects a wrong claim. The full-year RTS-GMLC_AC solve was reported as not completing. It hit openTEPES's 36,000 s limit at a barrier gap of 4.95e-09, essentially converged; Pyomo reports that as `aborted`, which was read as a failure. The log says "Time limit reached" four lines above it. The reason it needed the whole limit is memory, not difficulty. A one-month run was added to measure it: horizon rows factor barrier wall 1 day 54,855 - - 5.9 s 1 week 528,423 - - 97.5 s 1 month 2,266,671 2.4 GB 86.0 s 622.4 s 1 year 27,509,055 30.0 GB - 10 h cap Rows and factor memory scale linearly with the horizon. A 30 GB factor does not fit the 24 GB machine, so the year paged: Gurobi estimated 4 s per iteration and observed iterations took 500 to 1,400 s. The month is the largest horizon that fits, and its barrier takes 86 s of the 622 s wall; the rest is model build and result writing. With 40 GB or more the year should solve inside the existing cap. The case carries no candidates of any kind, so none of this involves investment decisions.
Adds a three-month run to the horizon table. Memory scales linearly with the horizon and the barrier scales linearly with rows, both measured: horizon rows factor barrier iters 1 month 2,266,671 2.4 GB 86.0 s 73 3 months 6,876,807 7.0 GB 255.3 s 79 1 year 27,509,055 30.0 GB - - Three months is 25.0% of the year's hours, 25.0% of its rows and 23.3% of its factor, so the 30 GB figure for the year is measured rather than extrapolated. The barrier does not degrade with size: 2.97x the time for 3.03x the rows, with the iteration count almost flat. Carried to the full year that is roughly 17 minutes of barrier with the factor resident. The 10 h solver cap is not what stops the year on this machine; 24 GB of RAM is. Wall time and the loose-branch count for the three-month run are still pending; it is in crossover.
openTEPES had one AC formulation, branch flow. Bose and Low prove the two SOC
relaxations give the same BOUND, and an earlier comment used that to justify not
writing the second one. The equivalence is a statement about the optimal value,
not about conditioning, solve time or behaviour inside branch and bound, so it
does not answer the question it was cited for.
IndACPowerFlow = 0 DC, as before
= 1 branch flow, as before
= 2 bus injection in W space, relaxed by a second-order cone
= 3 bus injection in rectangular coordinates, exact, needs a
non-linear solver
IndACCycle = 1 the loop condition around each independent cycle. Only
meaningful for 2: under branch flow the angle is a node
potential so the sum is identically zero, and the option is
refused there rather than silently ignored.
Measured on 24 hours of 9n_AC, with the same thermal limit, the same angle band
and the same tightened bounds:
branch flow 0.5462 MEUR reactive 3732 Mvarh Vmax 1.0024
bus inject. 0.5172 MEUR reactive 3967 Mvarh Vmax 1.0060
Of the remaining 5.3%, 4.5 points is the AC current penalty, which prices vCurr
and so exists only under branch flow. Remove it from both and the two agree to
about half a per cent. That is the equivalence measured rather than cited.
Three things were needed to get there, each worth knowing:
- The cone must be written in STANDARD form, ||(Wre, Wim, v)|| <= u with u and v
tied to the bus voltages by linear equalities. The natural rotated form leaves
an indefinite bilinear term whose convexity the solver has to infer, and the
inference is brittle: adding a constant to one side loses it and Gurobi then
reports the model as non-convex while QCPDual is set, which blocks its own
fallback.
- DualReductions must be off for modes 2 and 3, so apply_solver_options now takes
the model and can vary the preset by formulation.
- Gurobi's barrier is not reliable on this model: the same input solved twice and
failed once in three identical runs. ipopt solves it reproducibly to fifteen
figures, and PowerModels ships a separate conic model class for the same
reason. Mode 2 should be run on ipopt.
Also fixes a mistake this change introduced and then repeated. Guarding the whole
AC operation function on mode 1 left bus injection with no power balance; moving
the guard down left it with no reactive capability limit, no idle-unit gate, no
shunt definitions and no condenser gates, because those sit below it. Generators
then supplied nameplate Mvar untied to any active output and used eight times the
reactive power branch flow did. The branch flow relations are now wrapped rather
than returned out of, so nothing below can be skipped by accident.
Not working yet, and labelled as such in the module: mode 3 has never been run,
and the loop condition is linearised as sum(Wim/Vnom^2) = 0, which cuts the loop
mismatch about fourfold but moves the objective only in the eighth digit. The
tangent equality PowerModels uses in ACTPowerModel is the right form.
…ectangular IndACCycle now imposes Wim == tan(theta_i - theta_j) Wre, the form PowerModels uses in ACTPowerModel, instead of summing Wim over a constant around each cycle. vTheta is a node potential, so its differences sum to zero around any cycle by construction; tying the voltage product to it is what makes the recovered angles consistent. Measured on 24 hours of 9n_AC the worst loop mismatch goes from 9.3e-05 rad to 5.2e-18. The earlier version moved the objective only in the eighth digit because arg(W_ij) is Wim/Wre, not Wim over a fixed denominator. The cost barely moves, which is now the right answer rather than a symptom: once the reactive capability limits were restored the solution was already close to loop consistent, and the tangent makes it exact for nothing. Mode 3, rectangular, has been run for the first time. The four formulations line up as they should on the same 24 hours, relaxations below the exact optima: BIM W space 0.5171667 relaxation BIM W space + tangent 0.5171682 relaxation branch flow SOCP, no penalty 0.5180225 relaxation BIM rectangular 0.5302618 exact branch flow restored, no penalty 0.5305 exact The last two matter most. They reach the exact AC optimum by completely separate routes -- one imposes the exact current and angle relations on a fixed plan in branch flow space, the other solves the non-convex rectangular equations directly -- and agree to 0.04%. Different variables, different equations, different solver paths. That also closes the 10.5% gap this comparison started with: 5.6 points were reactive capability constraints that were not being built, 4.5 points the AC current penalty which exists only under branch flow, and about half a point genuine relaxation slack.
The two 7-day fixtures edited the bundled case in place and put it back in a finally. That works until the run does not reach the finally. A timeout, a Ctrl-C or a crash leaves the case truncated on disk, and it then gets worse: the next run reads the corrupted file as its own "original" backup and restores to that, so the damage becomes permanent and silent. It had already happened. sSEP was found with StageWeight 52 instead of 1 and a rewritten Duration table, left behind by a suite run killed on a timeout. Anyone who ran the suite and committed would have shipped it. Both fixtures now copy the case into tmp_path and edit the copy, so there is no restore step to miss. case_7d_binary already worked this way; the shared helper keeps the three from drifting apart. Checked by killing a run with SIGKILL part way through: no case file is modified, where before two were. Same 160 tests pass, and a completed suite now leaves the case folder clean.
IndACPowerFlow belongs in oT_Data_Option. A case that puts it in oT_Data_Parameter still built a full AC model, because the main read loop copies every Parameter column into the parameter dictionary. But the peek that runs before that loop only looked in Option, so it answered 0 and skipped the two AC-only tables, ReactiveDemand and BusShunt. The run then finished normally and reported a solved AC case that carried no reactive demand and no shunt devices. Only a warning marked it, so the result looked right and was not. The peek now checks both tables.
An existing shunt was wired in permanently: its injection is an equality, so it was in service every hour of the horizon. A case needing a capacitor bank out at light load could only delete the device, which also removed it at peak. A new Switchable column in oT_Data_BusShunt gives a device a state per hour. The candidate disjunction already had the right shape, with the investment decision as the state, so switchable and candidate devices now share one set of constraints and differ only in which variable plays that state. A device that is both is switched hourly and may close only once built. Two consequences worth naming. The reactive and active bounds now contain zero for a switchable device, the same trap candidates already avoided: a capacitor's range is otherwise strictly positive and it could never open. And the state is written to a new ShuntCommitment result, because a zero injection alone cannot distinguish an open bank from a closed one on a bus at zero volts. IndBinShuntSwitch picks a discrete state, the default and what a mechanically switched bank does, or a relaxed one that keeps an AC run continuous. The discrete form turns the AC second-order cone into a mixed integer cone, so the relaxed form is there for size. Devices are fixed unless the case says otherwise, so existing cases are unaffected. Measured on a 24 hour 9n_AC window with a 300 Mvar bank made existing: 0.5462 MEUR switchable against 0.6826 MEUR fixed in service, with the bank open in all 12 hours.
A shunt was a single device, in service or not. Real compensation is often a bank of identical units where the choice is how many are in, and that could not be said. A row with Units = N in oT_Data_BusShunt now becomes a bank of N identical devices, each carrying the full susceptance of one unit. The expansion happens when the table is read, so every set, variable, constraint and result downstream handles a bank with the machinery that already existed for single devices, and the number of units on is the position. This follows the VAR source model in Alvarez, Paredes and Rider, IET Generation, Transmission and Distribution 13(13), 2019, where a bus carries an integer count of sources of fixed susceptance rather than a continuous one. Two things here go further than that paper: the count is chosen per hour as well as per period, and the susceptance comes from the case instead of being fixed. Units of one bank are chained, in service and in the build decision. Three units of a four unit bank can otherwise be in service in four different ways at identical cost, and branch and bound walks all of them. The default is one unit, so existing cases are unaffected. Measured on a 24 hour 9n_AC window with a four unit bank, at three times the reactive demand: two units in one hour, three in another and four in the remaining ten, with no reactive slack.
Our AC gaps had nothing independent to be compared against. Every attempt to build a pure AC OPF case by hand produced an artefact instead: storage minimums left over from truncating a horizon, then capacitors sized at peak reactive demand and solved at an off-peak hour, which over-compensated and forced absorbing slack. pglib-opf removes the case construction from the question. It is CC-BY-4.0, it is the standard AC OPF benchmark, and it publishes a reference objective and a second-order cone gap for every case. The reader turns a MATPOWER case into a one hour openTEPES AC case with nothing on top of the network: no commitment, no reserves, no storage, no emissions, no investment. Anything else would enter the objective and make the comparison meaningless. It documents the two unit traps in the format, each worth a factor of 100. Bus Gs and Bs are MW and Mvar at one per unit voltage, not per unit, so the conversion divides by the base. Branch r, x and b are per unit already. On pglib_opf_case118_ieee, against a published AC OPF optimum of 97,214 dollars per hour: our exact rectangular model gives 97,100, and our second-order cone relaxation gives 96,624, a gap of 0.49% where pglib reports 0.91%. Ours is tighter because of the bound tightening pass. That places our formulations where the literature sits and shows the much larger gaps on RTS-GMLC_AC belong to that case, not to the model.
The case carried no shunt table at all, so the system had no reactive compensation of any kind. RTS-GMLC puts a 100 Mvar reactor on buses 106, 206 and 306, in RTS_Data/SourceData/bus.csv under "MVAR Shunt B", and zero everywhere else. That column is Mvar at one per unit voltage, so on this case's 100 MVA base each reactor is Bshb = -1.0. They are existing plant, not investment candidates and not switchable, exactly as the source has them. Adding them lowers the cost of the case from 61.65 to 60.00 MEUR. The direction is what the physics asks for: the 138 kV lines generate charging reactive, voltages sit against the 1.05 ceiling either way, and without the reactors that surplus had to be absorbed by generators at a cost. Every earlier measurement on this case, including the 61.65 figure, was therefore taken on a system missing its reactors. A test pins the three devices to the source values so the table cannot quietly drift or disappear again. RTS-GMLC_AC still has no shunt table and has the same gap.
The same gap the operational case had. RTS-GMLC puts a 100 Mvar reactor on buses 106, 206 and 306, and this case carried no shunt table, so it ran with no reactive compensation at all. Each reactor is Bshb = -1.0 on the case's 100 MVA base, existing plant and not switchable, exactly as the source data has them. The case is a full year, 8736 load levels, and is the one that previously took about ten hours and reached the memory ceiling, so it was not re-solved and there is no measured before and after cost for it. The operational case, which is small enough to solve twice, moved from 61.65 to 60.00 MEUR on the same change. What was checked here is that the three bound tightening tests that build this case still pass, so the reactors do not disturb the angle or voltage tightening. The reactor test now covers both cases instead of only the operational one.
Every RTS figure in sections 12 to 14 was measured on a system with no reactive compensation, because neither RTS case carried a shunt table. A new section 14.1 measures what the three reactors change on a 24 hour window: the cone goes from loose on 4 of 120 branches to tight on all 120, and the gap between the relaxed cost and the same plan under exact physics falls from 88% to 1.2%. Most of the looseness on RTS was the missing compensation, not the network. Section 12.2 read that looseness as a property of this network and now says so with a qualification. Sections 13 and 14 are marked withdrawn rather than restated. The RTS_SOC24 case they used was not kept, and a fresh window cut from RTS-GMLC_AC_Oper gives 8.4499 MEUR and 4 loose branches where they record 15.3720 and 9 of 120, on what should be the same hours of the same case. What differs has not been identified, so replacing their numbers with these would be wrong. The 168 hour figure is directly comparable and is restated, 61.65 to 60.00 MEUR. The changelog entries for this work are cut back to the house style, and the flag reading fix gains the entry it never had.
…ompare the four fairly The band was built for the W-space formulation only, so IndACPowerFlow = 3 ran with no angle limits at all. It is written in rectangular coordinates as a bilinear constraint, which the non-linear solver mode 3 already needs handles without trouble. The bug surfaced while re-deriving the formulation comparison, which had to be redone because every RTS measurement predated the reactors. A new section 14.2 gives the result on a 24 hour window with the reactors in place: all four formulations agree to within 0.07%, and branch flow and bus injection in W space differ by 0.003%, which is the Bose and Low equivalence as a measurement rather than a citation. Two corrections were needed before the table meant anything, and the second was the larger. AC_CURRENT_PENALTY is charged on vCurr, a variable only branch flow has. At its default the same window gives 8.06 MEUR for branch flow against 5.53 for rectangular, so branch flow appears 46% worse than an exact model it must bound from below. That is the penalty, not the formulation. Any comparison across formulations, or against another tool, has to zero it. The band is not binding on this window and fixing it changed nothing there. It also did not change the pglib case118 check, which still gives 97,100 dollars against a published 97,214. That 0.117% was expected to be the missing band, since pglib imposes 30 degrees, and is now ruled out by measurement. It remains unexplained.
A case can ask for one thing and get another, and until now the run finished and reported success either way. IndACPowerFlow written into oT_Data_Parameter rather than oT_Data_Option built a full AC model whose reactive demand and shunt tables were never read: 1438 Mvar of load silently became zero and only a warning marked it. The block now printed reads the values back off the BUILT model rather than off the case files, so a table that did not arrive shows up as a zero before the solve. Reactive demand and the shunt counts are reported as numbers for that reason, and an AC run with no reactive demand anywhere gets an explicit warning. It also notes that IndPTDF is a lossless representation, so a case asking for both PTDF and network losses is told that the losses are ignored. That is existing behaviour, not a change: the loss constraints have always been skipped when PTDF is on. It was simply never said.
The relaxation gap says whether the cone is tight. It does not say whether the operating point is physical, and those are different questions: a tight cone with a wrong branch equation passes the first and fails the second. That is exactly how the angle-relation sign error survived ten reviews and was caught only by an outside power flow. The check recomputes each branch flow from the bus VOLTAGES through the series relation and compares it with the flow the model reports. Deriving it from the flow variables would compare the flow equations with themselves and pass whatever they said. Until now this lived in prototypes/ and needed pandapower, which openTEPES does not ship, so a user could not run it at all. It agrees with that prototype exactly, 68.409915 MW from both on the relaxed 9n_AC solve. With IndACRestore = 1 the same case comes back at 0.00001 MW, which is the difference between a relaxed operating point and a physical one. A relaxed solve is not expected to sit on the series relation, so the residual is reported rather than judged. The new module also carries the susceptance matrices and the DC power transfer distribution factors. Those are not yet wired into the model: the factors are still read from oT_Data_VariablePTDF, and deriving them instead is a separate change. They are here because the residual check and the matrices share the branch enumeration, and assembling it twice would let the two copies disagree about which links are AC or what a tap means. The derivation is verified against the model's own DC formulation: on 9n it reproduces the angle constraint's flows to 0.000000 MW, and the flows it gives are independent of the reference node.
The AC feature had no coverage at all: none of its option flags, none of its input tables and none of its result files appeared in the published documentation. InputData gains the six AC flags in the options table, and sections for the two AC-only input files. The bus shunt section states the unit convention, which is the trap: Gshb and Bshb are per unit on SBase referred to the nominal voltage of the bus, so a device of -1.0 on a 100 MVA base is a 100 MVAr reactor. It also says that the injection follows the square of the voltage, so a bank delivers about 10% more than its nameplate at the upper voltage limit, and that each unit of a stepped bank carries the full susceptance of one unit rather than a share of the bank. OutputResults gains the twelve AC result files and the two diagnostics, with the distinction between them spelled out: the relaxation gap says whether the cone is tight, the residual says whether the operating point is physical, and neither substitutes for the other. The mathematical formulation page still has no AC content.
IndPTDF becomes an explicit three-valued option: 0 off, 1 reads the factors from oT_Data_VariablePTDF as before, and 2 computes them from the Reactance column. A case no longer has to produce them in another tool and paste in a branch-by-node table that openTEPES cannot check against its own network, and nothing can quietly disagree with the reactances beside it. The flag used to be implied by the presence of the table alone. That stays the default when a case says nothing, so existing cases are unaffected, but an explicit value now wins and a value contradicting the case is an error: asking to read a table that is not there, or to compute factors beside a table, both raise rather than picking one silently. The computed factors carry no load level index. The topology is fixed for a period, so an hourly index would store the same numbers once per hour: for a year of RTS-GMLC that is 8,736 x 120 x 73 entries of duplicated data. The read-from-table path keeps its hourly index, because a case may legitimately vary it. Mode 2 is refused when the case has candidate or switchable AC lines. The factors belong to one topology, and a decision that changes it would leave them stale in a way nothing detects: the flows stay plausible and stop being right. Generation candidates, storage candidates and candidate DC links all pass, because none of them enter the susceptance matrix. A case that does move its AC topology can still supply the factors itself. Verified against the formulation the factors stand in for: on 9n the computed ones give the same cost to the cent and the same branch flows to 0.000000 MW as the angle model. Two things had to be corrected on the way. pIndPTDF was declared within Binary, so a value of 2 would have been rejected outright. And four places read the PTDF table under a truthy test of the flag, which mode 2 reaches with no table to read; they now test for mode 1.
The page had no AC content at all: none of its parameters, variables or equations appeared in the published formulation. The new material follows the page's own organisation rather than standing apart from it. Parameters, variables and equations each gain an AC block in the section that already holds their kind, so a reader looking for all the parameters still finds them in one place. Writing it against the existing notation caught three collisions. The letter l is already the half ohmic losses, so the squared current is cu; f is already the active power flow, so the reactive one is q; and the voltage angle was already defined and is reused rather than redeclared. The reactance and the base power were also already there, and the series conductance and susceptance were used by the bus injection flows without being defined anywhere. A section appended at the end would have shipped all of that. Covers the reactive balance, the branch flow model, the two bus injection models, generator reactive capability, the shunt devices including the switchable and stepped ones, and the flow-based method. The prose is kept to one sentence per constraint, and only where the equation does not speak for itself: that the current definition is the sole relaxation, that the minus in the angle relation is deliberate, that the charging is left out of the bus injection flows to avoid counting it twice, and that the flow-based representation is lossless. Verified by building the documentation, which succeeds with no warnings.
…s together IndCycleFlow, IndSectorDecomposition, IndCompleteProblem and IndSequentialSolving were literals in openTEPES.py, so no case could select any of them. The model implements four stage-solving strategies, parallel, sequentially through an LP file, sequentially in memory and by sensitivity analysis, and all four were unreachable. IndSequentialSolving was also declared binary while its own code branches on those four values, so two of them were outside their own declared domain. They are now read from oT_Data_Option, with the values that used to be hard-coded as the defaults, so a case that says nothing behaves as before. Each of the four was run on 9n before being exposed: all give the same cost, and time Benders decomposition lands 8 EUR away, which is its convergence tolerance. Publishing a flag that leads to a broken path would be worse than leaving it hidden. The checks on incompatible options were made one at a time in five places, so a case with three clashing flags was told about one, fixed it, and was told about the next. They are collected into one pass and reported together, numbered. oT_Data_Option is now documented as the home of every indicator, with oT_Data_Parameter for the numeric scalars. Both are still read, because a flag reaching only one of them used to build a different model from the one the case asked for, but an indicator in the parameter file now says so. The run also states which problem and which stage strategy are in force. Named run-mode presets were considered and left out. With the flags exposed and the clashes reported together, a preset would be a second way of saying what the flags already say, and would raise a precedence question that does not exist today.
Our exact model reported 97,100 dollars per hour against a published 97,214, which is 0.117% below a figure nothing should beat. The cause is the form of the thermal limit, not a missing constraint. openTEPES writes the limit on the current, so with the cone the apparent power it admits is Smax times V_i over Vmin, which reaches Smax times Vmax over Vmin at the top of the voltage band. pglib imposes a flat cap on the apparent power. On case118 the band is 0.94 to 1.06, so the allowance is 1.1277, and the worst loaded branch of our solution sits at 112.8% of its rating. That is the whole of the difference. Scaling every rating by Vmin over Vmax, which makes the admitted apparent power never exceed the original, gives 97,431 and so brackets the published value from the other side. A constraint looser on one side and tighter on the other should do exactly that, which is why the bracket confirms the explanation rather than merely fitting it. Neither convention is wrong. A thermal limit is a heating limit and heating follows the current, so writing it on the current is the more physical of the two; a flat cap on apparent power is the more common. The consequence is only that an objective compared against pglib is not compared like for like unless the ratings are adjusted. The second-order cone gap is unaffected, since it is measured between our own relaxed and exact solves, which share the convention. The angle-difference band was the earlier candidate and had already been ruled out by measurement.
Sections 13 and 14 were withdrawn because every figure in them was measured on a system with no reactive compensation, and the case they used had not been kept. Rather than hunt for a case that no longer exists, both sections are measured again from scratch on a case whose definition is written down in the section, so this cannot happen a second time. Three things changed. The exact model does solve on RTS-GMLC. ipopt reports an optimal solution where the earlier attempt stopped at problem infeasible, and that attempt was made on a system with no reactors at all. It takes 1212 seconds against 10.5 for the cone, so the argument against using it as the working model is its cost, not any difficulty in solving it. The relaxation gap is 2.26% on the nine bus case and 0.003% on RTS. The small case is the loose one, which is the opposite of what section 12.2 concluded from measurements that included the current penalty. And the current penalty is what makes the cone look tight. Charging for the current pins the relaxation to the cone boundary, because the cheapest current consistent with the flows is the one on it. With the penalty on, the relaxed and exact costs of the nine bus case agree to six figures and the cone reports itself tight; with it off, the same case has a 2.26% gap and a cone loose on two thirds of its branches. An earlier reading of that agreement as "where the cone is tight the relaxed answer is the AC answer" was wrong: the agreement was real, but what it showed was the penalty doing its work. Section 14 now carries three costs rather than two, because the relaxation understating its own plan by 0.55%, the relaxation gap of 0.003%, and the relaxed plan being 0.55% worse than the best available are three different statements. Two notes made stale by this rewrite, in the section 12 header and in 14.1, are corrected.
Section 12 was the last part of these notes still carrying figures taken on a system with no reactive compensation. Two of the new numbers check the change rather than merely report it. The DC model is unchanged to the digit, which is what should happen because a shunt has no place in a DC model and never reaches it. The AC model grew by exactly 504 rows and 504 columns over the week, which is 3 reactors over 168 hours. Neither was arranged. The relaxation tightness is the substantive change. The reactors cut the worst cone gap from 0.901 to 0.120, a factor of seven, and left more branches marginally loose, 31 against 21. Severity down and spread up: compensating the system stops any one branch being badly wrong and leaves a larger number slightly wrong. The earlier headline, that the loose cone on RTS is real, still stands, but most of its magnitude was missing reactive support. A third row is added with the current penalty off, where two thirds of the branches are loose. Section 13.1 explains why: the penalty charges the current and so pins the relaxation to the cone boundary, which means the first two rows measure the penalty as much as the network. The horizon table keeps its three month and full year rows for their sizes and marks them as not re-solved, since they need 7 GB and 30 GB. Each row now names the case it came from. The old 24 hour row recorded 54,855 rows where the same horizon now gives 75,087, because that row came from the RTS_SOC24 case that was not kept, which is the same case behind the figures withdrawn from sections 13 and 14. Wall times are marked indicative: the machine was doing other work.
Commit 341f223, titled "modify comments", also removed eighteen lines of working code from openTEPES_ProblemSolvingResolve: the module-level import of enabled and resolve_persistent, and the whole opt-in block in resolve that hands the re-solves to the persistent solver. The effect was that --warm-resolve, added on purpose in e238c01 as a persistent solver for Mode C sweeps, has done nothing since 4 August. Switching it on fell through to the standard path with no message, so the only symptom was that the re-solves stayed slow. The warm sweep module itself was untouched; only the wiring to it was cut. test_mode_c_resolve_warm_resolve_parity was the only thing detecting this and had been failing ever since. The test was right and is not changed. The restored block is the text from 341f223^, not a reconstruction, and it goes back in both branches of the try/except import so the direct-run path works as well.
The penalty prices the branch current so the relaxation cannot buy voltage
with current that is not there. It is a numerical device, not money, and it
was inside the reported total: on a 168 hour RTS-GMLC window it came to
14.43 MEUR of a 60.00 MEUR figure, a quarter of it. The code comment put it
at about 2.4%, which is what it measures on the nine bus case; RTS is an
order of magnitude worse. It also reached the locational prices through the
eBalanceElec duals.
It now has its own variable, which enters the OBJECTIVE but not
vTotalSCost:
objective = vTotalSCost + sum(scenario factor * vTotalNPenalty)
vTotalSCost = the costs, with no penalty in them
The solver is steered exactly as before and every writer, summary and test
that reads vTotalSCost gets money rather than money plus a device, without
one of them having to change.
The split reconciles: the same case now reports 45.5633 MEUR with a
separate 14.4333, and those sum to the 59.9967 it reported before. The
objective value is identical, so the optimum has not moved, only the way it
is broken down.
The economic results no longer subtract the penalty back out of the network
operation cost, because vTotalNCost never carries it now, and the row is
renamed so that nobody adds it into the total.
This does change a user-visible number: an AC case reports a lower total
than it did, and the difference is the device that was never a cost.
erikfilias
marked this pull request as draft
August 21, 2026 18:10
The unit CI job took between twenty-four and fifty-two minutes without solving anything. Most of that was five builds of RTS-GMLC_AC, the full year case, at 70.5 s each against 2.0 s for the operational case. Four of those five were the bound tightening tests, and they do not need a year. The two cases carry the SAME network: 73 nodes, 120 branches, and tightened angle and voltage bounds that compare bit-identical, because the operational case is a week of the same system. Tightening is a property of the network and not of the horizon. That is recorded in a comment above the section, so the full year case is not restored later by someone assuming coverage was lost. The fifth build stays. The reactor test is parametrised over both RTS cases precisely to check that each carries the three reactors, so pointing it at one of them would delete half of what it asserts. Measured on this machine: the unit selection drops from most of half an hour to 5m47s, the AC suite from 10m56s to 6m02s, and the whole suite to 9m20s. Unrelated, and not fixed here: one run of the full suite failed test_every_formulation_solves_and_writes_its_results[2-0-0-ipopt], and two later runs did not reproduce it. That is the bus injection formulation in W space on a non-convex solve, whose fragility is already recorded in the module header. It is a flaky test and will show up in CI eventually.
…e prose Two gaps found by checking the pages against the changelog. The parameter table listed thirteen entries and none of the AC ones. VMin, VNom, VMax, CapacitivePF, InductivePF, ConverterPF and EpsilonCurrent are all read from oT_Data_Parameter and used by the AC code, and none appeared. ConverterPF is not in any bundled case either, so a case turning on IndACConverter had neither an entry to read nor an example to copy. The HVDC converter model had no equations. IndACConverter was listed among the options, but the model itself was not written down: the line-commutated converter drawing reactive power at both terminals through the split flow, and the voltage-source converter bounded by the converter rating and released only in service. Both are added, with the two variables they use. The narrative is also made plainer. Several sentences stated a view rather than a fact: that a relaxation should never exceed the exact optimum, that an AC case with no reactive demand is almost always a mistake, which file to read first, what a wrong branch equation looks like. They are replaced with what the files contain and what the numbers were. One of them also said "two files" above a table of three. The documentation builds with no warnings.
Active and reactive power were bounded separately, so a station could hold P at the link rating and Q at tan(acos(pf)) times it at the same moment, delivering NTC/pf of apparent power. At the default power factor of 0.85 that is 17.6% more than the converter has; at 0.95 it is 5.3%. The rating is a rating on the apparent power, so it now bounds the apparent power. It is written as a ring of twelve tangent lines rather than as the disc, for two reasons found by testing rather than by argument. A quadratic constraint would put a cone into IndACModelType = 1, whose purpose is to stay a mixed-integer linear problem an LP/MIP solver can take, and that mode is the only AC variant the project's CI can run at all. And on a tightly rated link the disc made the barrier stop with "numerical trouble" on a case that solved without it. The cuts circumscribe the disc, so the bound is loose by 1/cos(pi/12), which is 3.5%. Measured: a voltage-source converter settles at 103.53% of its rating and a line-commutated one at 100.00%, the latter exact because its reactive power is a fixed ratio of its active power and the limit collapses to a linear bound on the active flow. Trading 17.6% of free capacity for 3.5% of looseness, and keeping the MILP path, seemed the better side of that trade. CONV_CUTS is one line if a case needs it tighter. Also adds ConverterPF to 9n_AC, which had no case carrying it, and a test that the two converter models move the system cost in opposite directions: a line-commutated station draws reactive power and costs the system, a voltage-source station supplies it and relieves the system. The existing converter tests check the mechanism and would not catch a reversed sign.
The price on the branch current is in the objective, so it reaches the duals of the nodal balance. Section 17 measures it on a 24 hour RTS window. The distortion is not a level shift: removing the mean movement leaves the worst deviation where it was. The prices move relative to one another, which is what a comparison of nodal against zonal pricing reads. At the value the code carried before the price became case data the worst nodal price moved by 140 EUR/MWh against a spread of 200; at the value the RTS cases now carry it is 1.53. Section 17.1 records the attempt to take prices from the restoration pass instead. The mechanism works, the values did not check out, and the change was reverted. Written down so the next attempt starts from the open question rather than from the experiment.
doc/ is for the published documentation. The working notes from building the AC optimal power flow were sitting alongside it and are not part of what the model documents, so they move to the project repository that the work serves. Nothing users read changes. The feature is documented in doc/md, which is untouched.
The diagram was stale in three layers. Layer 4 said six formulation files and there are eight. Layer 6 said twelve results modules and there are thirteen. Layer 3 did not name the modules the AC path adds. Layer 4 gains an ELECTRICITY NETWORK bracket under the electricity box holding the three interchangeable network models: DC, branch flow and bus injection. The box it replaces was the planned one for selectable dc_opf / ac_opf builders, so a planned item becomes an implemented one. The AC files sit under electricity rather than beside Hydro, Hydrogen and Heat, because they are a choice of network model and not a new energy carrier. The PNG is re-rendered from the SVG at the width doc/img/README.md documents.
vTheta appears in no constraint when the loop condition is off, so the solver sets some nodes and leaves others unset, and which ones varies between runs and between solvers. The guard that decides whether a nodal voltage phasor can be formed returned on the first node it found. Whenever that node happened to carry a value it reported the angles available, and the residual check then built a phasor from None at a later node and raised TypeError. It now tests every node. This is what made IndACPowerFlow = 2 with IndACCycle = 0 fail once in every few runs of the test suite. The new test fails against the old guard and passes against this one, so it pins the bug rather than the symptom.
A converter station is not free to pass power through, and openTEPES modelled an HVDC link as active power transfer with a line loss factor and nothing at the terminals. For a study of nodal against zonal pricing that matters: what a link costs to cross is part of the price spread across it. Each terminal carries a station and each station is charged separately, which is the convention the existing DC line loss factor already uses. The no-load part is drawn while the link is in service, as a fraction of the link rating. The marginal part is drawn on the power the station carries, either direction. Set them with ConverterNoLoadLoss and ConverterMarginalLoss; both default to zero, so a case that does not ask for them gets the results it got before. The marginal part needs the size of the link flow, and getting that exactly needs the flow-direction binary the line-commutated model already carries. A loss therefore brings that binary in under the voltage-source model too. Without it the model can inflate both halves of the flow and discard surplus energy into a loss that does not exist, which is the same failure the reactive draw had. Reported per link in oT_Result_NetworkConverterLosses.
The page background stopped 45 px short of the canvas, and the arrow from resolve.py back to SettingUpVariables.py was routed through that unpainted strip, so it appeared to run off the edge. The arrow now stays inside the page. Its rotated label is gone: it repeated the sweep-modes panel word for word, and removing it freed the space the arrow needed. The footer ran off both edges and lost its opening words. It is shorter now and fits. Claims about the design gave way to descriptions of it. "single source of truth" is "column and type specs", "drop-in backend" is "same InputSource interface", and "cost: solve only (cheapest)" drops the judgement. The note about splitting the other sectors moved back into the planned box, where it belongs, having been left under an implemented group by the previous change. The planned boxes were checked against the code. All five are still planned.
Each bend was one quadratic curve spanning its whole segment, with the control point at the corner, so the curve radius followed the segment length rather than being chosen. One bend swept over 435 px and the other over 27 px, and they did not read as a pair. The path now has straight runs with a quarter turn of equal radius at each corner.
HiGHS cannot express a nonlinear constraint, so the second-order cone, the exact non-linear model and the AC restoration pass were skipped on every runner. The cone is the default, which meant the formulation carrying every validation number was never solved in CI. A conic solver is not needed for this. The relaxation is convex, so ipopt reaching a local optimum reaches the global one; on the 9n case gurobi and ipopt agree to 3e-07 relative. The formulation matrix gains a branch-flow cone row under ipopt, because the existing cone row asks for gurobi and there is no licence for it on a runner. ipopt joins the conda transaction that already installs the test tools, so it costs one solve rather than two. The tests it unlocks run in under ten seconds. Linux only, because ipopt's convergence depends on how MUMPS was built and a platform-dependent solver makes for flaky tests. A new step prints the solvers each runner has, so the next gap is visible rather than silent. This tests that those paths build and solve. It does not stand in for a conic solver certifying the bound.
It asked for gurobi, and the bundled licence is size-limited, so it was skipped everywhere and the restoration pass had no CI coverage at all, even after ipopt arrived on the Linux job. The outer solve is the second-order cone, which is convex, so an NLP solver reaches the same optimum. The pass itself already ran on ipopt.
The AC work introduced seven new files. Each one now sits in the module that already owned that concern, so the package has the same file list as master: reading and bound tightening -> openTEPES_InputData.py network matrices and AC set-up -> openTEPES_DataConfiguration.py AC variables -> openTEPES_SettingUpVariables.py branch flow, bus injection, converters, restoration -> openTEPES_ModelFormulationElectricity.py AC results -> openTEPES_OutputResultsNetwork.py Function names and their callers are unchanged; only the file they live in moved. Top-level names were checked for collisions before the two formulation modules were merged into one, and there were none. The architecture diagram advertised the seven modules, so it is corrected: the formulation layer is six files again, and the three network models are shown inside Electricity.py. The middle box is labelled BF rather than AC, which was wrong beside BIM because bus injection is an AC model too. The functions keep their AC names, correctly: they run for every AC mode and only their interior is gated on branch flow.
# Conflicts: # CHANGELOG.md
NetworkReactiveNotServed spelled out what every other not-served result abbreviates: NetworkENS, NetworkPNS, NetworkHNS. It is NetworkQNS now, which is also what the variable behind it has always been called, vQNSPos and vQNSNeg. NetworkUtilizationAC put the carrier after the word, where NetworkElecUtilization and NetworkHeatUtilization put it before. It is NetworkElecUtilizationAC. Neither name has been released, so nothing downstream depends on them. This also restores three changelog entries that were lost when an earlier commit was split: the relocation, the network matrices wording, and the corrected description of the diagram.
The AC section of the mathematical formulation had no figures, and it is the
most notation-dense part of the documentation. Three single-line diagrams now
sit next to the equations they describe:
ac_branch_model one branch: the tap, the series impedance and half
the charging susceptance at each end, with the
sending and receiving flows at their own ends
ac_hvdc_converter the link under both converter models, with the
reactive power and the station losses at each
terminal
ac_converter_capability the capability disc of one terminal and the twelve
tangent lines that stand in for it
Symbols follow the notation table. Parameters are drawn in blue and variables in
red, so a reader can tell what the model decides from what the case file fixes.
The arrangement matches the diagrams already in the repository: hand-drawn SVG
as the source, a PNG rendered from it, and the documentation embeds the PNG.
doc/img/README.md says how to regenerate them.
The captions leaned on capitals for emphasis and on connectives that argued a point rather than stating one. DRAW, EACH, SUPPLIES, ABSORBS and OPPOSITE are lower case now, and clauses of the form "which is why" are gone. Four captions became one sentence each. The angle relation stands on its own without a gloss naming what it ties together, and the sentence explaining why a loss factor is unnecessary is replaced by the fact it rested on: the sending and receiving flows are separate variables and their sum is the loss.
The F-bar label sat on the dashed radius, so the line ran through the text. The radius now points into the upper-left quadrant, which is empty, and the label sits beside it rather than on it.
The vertical axis ran to y=465 and the first caption sits at y=452, so the line crossed the text. The polygon reaches y=421, so the axis stops at 434 and still extends past everything it has to.
The separator in the right-hand converters ran the same way as the left, from bottom-left to top-right, but their AC and DC marks sit in the opposite corners because the DC side faces the link. Both marks therefore landed on the line. The separator is mirrored in both rows, so each mark sits in a corner the line avoids. The two stations of a link now mirror each other, with AC facing its own node and DC facing the link.
The README described the operational model as network-constrained unit commitment via DC power flow, and the ohmic losses as proportional to the line flow. Both read as unconditional and are no longer so. An AC power flow is named as an alternative that is off by default, the losses follow the exact relation under it, and the result topics list the voltage magnitudes, reactive flows, shunt injections and reactive-power marginal that come with it. Four sentences changed; nothing else moved.
One bug, four likely, four readability. Each was checked against the tree by
content, because the line numbers in the report predate the module relocation.
B1 The stage-solve path builds its own call list rather than walking
FORMULATION_REGISTRY, and named two of the three AC blocks. Bus injection
was missing, and NetworkACCurrent returns without doing anything outside
mode 1, so on that path modes 2 and 3 had their balances with nothing tying
the flows to the voltages: a free transport network reporting a lower cost
and no message. A test now compares the registry against that path's source,
which catches the drift without a Benders solve.
L1 The stage objective left out the price on the branch current, so vCurr was
unpriced there and the relaxation could buy voltage with current that does
not exist.
L2 eBIMSLimit gates the cone on the line state, which is a product of two
variables when that state is free. Candidate and switchable AC lines are now
refused under bus injection, as they already are under IndPTDF = 2. The
split-cone alternative that would lift this is recorded in the comment.
L3 The angle writer read vTheta without the guard the diagnostic in the same
module already uses, so it could build a phasor from None.
L4 A condenser read InvestmentUp = 0 as "not buildable" while generators,
network, H2 pipes, heat pipes and bus shunts all read it as "no limit". The
condenser now follows the same convention. No shipped case has a condenser,
so this has no test.
R1 A comment described a version of eBIMSLimit that is not the code, and said
the opposite of what the constraint does.
R2 Two package imports sat outside the try/except the rest of the package uses.
R3 bsh was computed once per row of four constraints and never read.
R4 The restoration assumes the one active objective is the system one, which
the Benders path need not leave behind. Noted where it matters.
vTotalNPenalty was declared for every case and added to the objective unconditionally, but the constraint that defines it is skipped when AC is off. A DC model therefore carried one unconstrained column per load level, priced in the objective. Being non-negative and minimised, they settle at zero, so the reported cost was always right; they are dead columns all the same, and they change what the solver presolves. Both the variable and the objective term now exist only under IndACPowerFlow = 1. On 9n_H2 the variable count now matches the base tree exactly. The cost summary had also gained two AC rows on every case, both zero where there is no AC power flow, which changed the shape of a file every case writes. They are written only when AC is on. Neither changed the objective: 9n_H2 reports 168.3419672985 before and after, and on the base tree. That case reaches a different dispatch at the same cost from one run to the next on the base tree too, so a file-by-file comparison of its results says more about degeneracy than about any change.
The review read the arithmetic of TightenACBounds and found it consistent, but
could not reproduce it numerically, and named the three-node triangle as the
test that would close it.
Three nodes in a ring, the first of them the reference, with a 1.03 tap on one
branch so its effect is visible against two untapped twins:
the closed form the bound is asin(Smax z / (Vmin tau Vmin)), checked
against arithmetic done in the test, and the tapped
branch comes out tighter than its twins, which is what
puts the tap on the sending voltage and nowhere else
the declared band a band narrower than the implied one survives untouched,
and a one-sided band is not made symmetric
candidate branches with every branch a candidate, no bus is tightened,
because a released drop equation implies nothing
Each was checked against a mutation: dropping the tap from the divisor fails the
first, and letting candidates propagate fails the third.
The condenser is modelled — recognised from the data, gated on its investment decision, its reactive output bounded by what the case declares, its cost inside the electricity fixed cost — but no shipped case has one, so none of that path ran. Two tests build one at Node_5 of 9n_AC. The first checks an existing condenser is recognised and absorbs reactive power within its declared band. The second builds a candidate whose InvestmentUp is 0 and checks it comes out buildable, which is what that column means for generators, network, H2 pipes, heat pipes and bus shunts. It fails against the reading the condenser had before. An earlier commit said this could not be tested because no case has a condenser. That was the wrong conclusion: no case has one yet, and building one is thirty lines. The trap is that a unit must be declared in oT_Dict_Generation, and its technology in oT_Dict_Technology, not only added to the data table, or it never enters the generator set.
The condenser tests went in green and skipped on every runner: they asked for gurobi, whose bundled licence is size limited, so the path they were written for gained no CI coverage at all. Checking the runner log rather than the tally is what showed it. Four move to ipopt, which a runner has since the Linux solve job gained it: both condenser tests, the voltage-source converter, and the angle guard. None needs a binary, and the case each solves is the cone, which is convex. Five converter tests still ask for gurobi and still skip there — the LCC draw, the station losses, the apparent power rating, the two-model cost comparison and the unbuilt HVDC candidate. Each needs the flow-direction binary, and ipopt relaxes it: the LCC test fails outright when switched. Covering those needs either a mixed-integer variant of each on the piecewise model type with HiGHS, or a solver licence in CI.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Adds an AC optimal power flow to openTEPES, off by default.
Of 348k added lines, 341k are case data: 182 CSV files, mostly 8736-hour series for the new RTS-GMLC AC cases. The code is +4,459 lines and the documentation +619. No new modules: each concern sits in the file that already owned it.
Options and formulations
IndACPowerFlowselects the network model: 0 DC, the default and unchanged; 1 branch flow; 2 bus injection in W space; 3 bus injection in rectangular coordinates.IndACModelTypechooses how the branch current is written, as a second-order cone, a piecewise staircase, or the exact non-linear equation.IndACRestorere-solves the network at the exact equations to recover a physical operating point.It also adds generator reactive capability and synchronous condensers, bus shunts with an hourly on/off state (
Switchable) and stepped banks (Units), HVDC converter models with a power factor, an apparent power limit and station losses (IndACConverter), and voltage and angle bound tightening.Separately from the AC work,
IndPTDFbecomes three-valued and can compute the distribution factors from the reactances instead of reading them from a table.Validation
Against pglib-opf case118, whose optimum and SOC gap are published:
The −0.117% is a convention difference. The thermal limit here is written on the current, so it admits
Smax·V/Vminwhere pglib caps apparent power flat; scaling the ratings byVmin/Vmaxbrackets the published value from the other side at +0.223%.On a 24-hour RTS-GMLC window the four formulations agree to within 0.07%, and branch flow and bus injection differ by 0.003%.
The flows are checked against pandapower to 2e-09 p.u., and by a residual check in the model that recomputes each branch flow from the bus voltages: 68 MW off on a relaxed solve, 0.00001 MW after restoration.
Changes to existing behaviour
vTotalSCost. Its value is case data (EpsilonCurrent), because it is not scale-invariant.bus.csv. Adding them moves the operational case from 61.65 to 60.00 MEUR and cuts the worst relaxation gap sevenfold.Tests and continuous integration
Full suite: 189 passed.
CI runs HiGHS, which cannot express a nonlinear constraint, so the Linux solve job also installs ipopt. That covers the second-order cone, the exact non-linear model, both bus injection formulations and the restoration pass. The cone is convex, so ipopt reaches the same optimum a conic solver would; on the 9n case the two agree to 3e-07. It is Linux only, because ipopt's convergence depends on how MUMPS was built.
This shows those paths build and solve. It does not stand in for a conic solver certifying the bound, which is done locally and by the pglib check above.
A week-long AC case is not viable in CI: HiGHS takes 9,189 s on 168 hours of RTS-GMLC against a thirty-minute per-test timeout.
test_run.pysolves 9n_AC over 24 hours instead, which exercises the same writers in about three seconds.