Cannot reproduce IEEE-9 loadflow results

[fhjgch]

Background

For an experiment, we want to reimplement this grid, which follows the IEEE9 node example:

The grid specifications are given as follows:

According to this text book, the resulting power flow solution should be:

Problem Description

When we try to reimplement this grid in GridCal, using the attached code, we get this output:

# ieee9 grid
grid9 = gce.MultiCircuit('IEEE-9', Sbase=100)

# Add the buses and the generators and loads attached
bus1 = gce.Bus('Bus 1', Vnom=17.16, is_slack=True)#
grid9.add_bus(bus1)
gen1 = gce.Generator('Slack Generator', P=0.0, vset=0.0746)#14.145#1.0
grid9.add_generator(bus1, gen1)
bus2 = gce.Bus('Bus 2', Vnom=18.45)#18.45
#grid9.add_load(bus2, gce.Load('Gen 2', P=-163))
gen2 = gce.Generator('Gen2', P=163, Sbase=100, vset=0.08022)#14.145
grid9.add_bus(bus2)
grid9.add_generator(bus2, gen2)
bus3 = gce.Bus('Bus 3', Vnom=14.145)#
grid9.add_bus(bus3)
gen3 = gce.Generator('Gen3', P=85,  vset=0.0615)#
#grid9.add_load(bus3, gce.Load('Gen 3', P=-85))
grid9.add_generator(bus3, gen3)
bus4 = gce.Bus('Bus 4', Vnom=230)
grid9.add_bus(bus4)
bus5 = gce.Bus('Bus 5', Vnom=230)
grid9.add_bus(bus5)
grid9.add_load(bus5, gce.Load('Load 1', P=125, Q=50))
bus6 = gce.Bus('Bus 6', Vnom=230)
grid9.add_bus(bus6)
grid9.add_load(bus6, gce.Load('Load 2', P=90, Q=30))
bus7 = gce.Bus('Bus 7', Vnom=230)
grid9.add_bus(bus7)
bus8 = gce.Bus('Bus 8', Vnom=230)
grid9.add_bus(bus8)
grid9.add_load(bus8, gce.Load('Load 3', P=100, Q=35))
bus9 = gce.Bus('Bus 9', Vnom=230)
grid9.add_bus(bus9)

# add Lines connecting the buses
tr1 = gce.Transformer2W(bus4,
                                         bus1,
                                         'T1',
                                         HV=230,
                                         LV=16.5,#,
                                         short_circuit_voltage=14.3,
                                         r=0,
                                         x=0.2889,
                                         rate=247.5,
                                         tap_phase=0)#2.618
tr2 = gce.Transformer2W(bus7,
                                         bus2,
                                         'T2',
                                         HV=230,
                                         LV=18,
                                         short_circuit_voltage=12,
                                         r=0,
                                         x=0.3125,
                                         rate=192,
                                         tap_phase=0)#2.618
tr3 = gce.Transformer2W(bus9,
                                         bus3,
                                         'T3',
                                         HV=230,
                                         LV=13.8,
                                         short_circuit_voltage=7.5,
                                         r=0,
                                         x=0.293,
                                         rate=128,
                                         tap_phase=0)#2.618
grid9.add_transformer2w(tr1)#0.5236 => 30°#2.618
grid9.add_transformer2w(tr2)#2.618#2.618#0.5236
grid9.add_transformer2w(tr3)#0.5236
grid9.add_line(gce.Line(bus_from=bus4,
                        bus_to=bus5,
                        name='line 4-5',
                        r=0.050095,
                        x=0.425331,
                        b=0.035232,
                        contingency_enabled=False))
grid9.add_line(gce.Line(bus_from=bus4,
                        bus_to=bus6,
                        name='line 4-6',
                        r=0.085066,
                        x=0.460302,
                        b=0.031576,
                        contingency_enabled=False))
grid9.add_line(gce.Line(bus_from=bus5,
                        bus_to=bus7,
                        name='line 5-7',
                        r=0.159735,
                        x=0.805293,
                        b=0.061158,
                        contingency_enabled=False))
grid9.add_line(gce.Line(bus_from=bus7,
                        bus_to=bus8,
                        name='line 7-8',
                        r=0.194707,
                        x=0.849716,
                        b=0.071462,
                        contingency_enabled=False))
grid9.add_line(gce.Line(bus_from=bus6,
                        bus_to=bus9,
                        name='line 6-9',
                        r=0.042533,
                        x=0.360113,
                        b=0.028917,
                        contingency_enabled=False))
grid9.add_line(gce.Line(bus_from=bus8,
                        bus_to=bus9,
                        name='line 8-9',
                        r=0.059546,
                        x=0.503781,
                        b=0.041880,
                        contingency_enabled=False))
## Powerflow
options = gce.PowerFlowOptions(gce.SolverType.NR,
                               retry_with_other_methods=False,
                               tolerance=1e-4,
                               control_q=False,
                               control_taps_phase=False,
                               control_taps_modules=False,
                               apply_temperature_correction=False,
                               verbose=False)


drv = gce.InputsAnalysisDriver(grid=grid9)
mdl = drv.results.mdl(gce.ResultTypes.AreaAnalysis)
df = mdl.to_df()
print(df)



power_flow = gce.PowerFlowDriver(grid9, options)
power_flow.run()
print(grid9.name)
print('Converged:', power_flow.results.converged, 'error:', power_flow.results.error)
resultsDF = power_flow.results.get_bus_df()
print(resultsDF)

We see that the voltages at the buses are far from the expected result. What are we missing in our implementation in order to reproduce the results from the book?
Thank you very much in advance!

[SanPen]

Hi!

The best way I'd model that example is:

# https://github.com/SanPen/GridCal/issues/354
import numpy as np
import GridCalEngine.api as gce

# ieee9 grid
grid9 = gce.MultiCircuit('IEEE-9', Sbase=100, fbase=50)

# Add the buses and the generators and loads attached
bus1 = gce.Bus(name='Bus 1', Vnom=17.16, is_slack=True)
bus2 = gce.Bus(name='Bus 2', Vnom=18.45)
bus3 = gce.Bus(name='Bus 3', Vnom=14.145)
bus4 = gce.Bus(name='Bus 4', Vnom=230)
bus5 = gce.Bus(name='Bus 5', Vnom=230)
bus6 = gce.Bus(name='Bus 6', Vnom=230)
bus7 = gce.Bus(name='Bus 7', Vnom=230)
bus8 = gce.Bus(name='Bus 8', Vnom=230)
bus9 = gce.Bus(name='Bus 9', Vnom=230)

grid9.add_bus(bus1)
grid9.add_bus(bus2)
grid9.add_bus(bus3)
grid9.add_bus(bus4)
grid9.add_bus(bus5)
grid9.add_bus(bus6)
grid9.add_bus(bus7)
grid9.add_bus(bus8)
grid9.add_bus(bus9)

grid9.add_generator(bus1, gce.Generator(name='Slack Generator', P=0.0, vset=1.0))
grid9.add_generator(bus2, gce.Generator(name='Gen2', P=163, Sbase=100, vset=1.0))
grid9.add_generator(bus3, gce.Generator(name='Gen3', P=85, vset=1.0))

grid9.add_load(bus5, gce.Load(name='Load 1', P=125, Q=50))
grid9.add_load(bus6, gce.Load(name='Load 2', P=90, Q=30))
grid9.add_load(bus8, gce.Load(name='Load 3', P=100, Q=35))

# add Lines connecting the buses
tr1 = gce.Transformer2W(bus_from=bus4,bus_to= bus1, name='T1', HV=230, LV=16.5, nominal_power=247.5, rate=247.5, tap_phase=150*np.pi/180)
tr1.fill_design_properties(Pcu=0.0, Pfe=0.0, I0=0.0, Vsc=14.3, Sbase=grid9.Sbase)

tr2 = gce.Transformer2W(bus_from=bus7, bus_to=bus2, name='T2', HV=230, LV=18, nominal_power=192, rate=192, tap_phase=150*np.pi/180)
tr2.fill_design_properties(Pcu=0.0, Pfe=0.0, I0=0.0, Vsc=12.0, Sbase=grid9.Sbase)

tr3 = gce.Transformer2W(bus_from=bus9, bus_to=bus3, name='T3', HV=230, LV=13.8, nominal_power=128, rate=128, tap_phase=150*np.pi/180)
tr3.fill_design_properties(Pcu=0.0, Pfe=0.0, I0=0.0, Vsc=7.5, Sbase=grid9.Sbase)

grid9.add_transformer2w(tr1)  # 0.5236 => 30°#2.618
grid9.add_transformer2w(tr2)  # 2.618#2.618#0.5236
grid9.add_transformer2w(tr3)  # 0.5236

l1 = gce.Line(bus_from=bus4, bus_to=bus5, name='line 4-5')
l1.fill_design_properties(r_ohm=5.3, x_ohm=45.0, c_nf=1060, length=1.0, Imax=1.0, freq=grid9.fBase, Sbase=grid9.Sbase)

l2 = gce.Line(bus_from=bus4, bus_to=bus6, name='line 4-6')
l2.fill_design_properties(r_ohm=9.0, x_ohm=48.7, c_nf=950, length=1.0, Imax=1.0, freq=grid9.fBase, Sbase=grid9.Sbase)

l3 = gce.Line(bus_from=bus5, bus_to=bus7, name='line 5-7')
l3.fill_design_properties(r_ohm=16.9, x_ohm=85.2, c_nf=1840, length=1.0, Imax=1.0, freq=grid9.fBase, Sbase=grid9.Sbase)

l4 = gce.Line(bus_from=bus6, bus_to=bus9, name='line 6-9')
l4.fill_design_properties(r_ohm=20.6, x_ohm=89.9, c_nf=2150, length=1.0, Imax=1.0, freq=grid9.fBase, Sbase=grid9.Sbase)

l5 = gce.Line(bus_from=bus7, bus_to=bus8, name='line 7-8')
l5.fill_design_properties(r_ohm=4.5, x_ohm=38.1, c_nf=870, length=1.0, Imax=1.0, freq=grid9.fBase, Sbase=grid9.Sbase)

l6 = gce.Line(bus_from=bus8, bus_to=bus9, name='line 8-9')
l6.fill_design_properties(r_ohm=6.3, x_ohm=53.3, c_nf=1260, length=1.0, Imax=1.0, freq=grid9.fBase, Sbase=grid9.Sbase)

grid9.add_line(l1)
grid9.add_line(l2)
grid9.add_line(l3)
grid9.add_line(l4)
grid9.add_line(l5)
grid9.add_line(l6)

gce.save_file(grid9, "IEEE9.gridcal")

# Power flow
options = gce.PowerFlowOptions(gce.SolverType.HELM,
                               retry_with_other_methods=False,
                               tolerance=1e-6,
                               control_q=False,
                               control_taps_phase=False,
                               control_taps_modules=False,
                               apply_temperature_correction=False,
                               verbose=False)

drv = gce.InputsAnalysisDriver(grid=grid9)
mdl = drv.results.mdl(gce.ResultTypes.AreaAnalysis)
df = mdl.to_df()

print(df)

power_flow = gce.PowerFlowDriver(grid9, options)
power_flow.run()
print(grid9.name)
print('Converged:', power_flow.results.converged, 'error:', power_flow.results.error)
resultsDF = power_flow.results.get_bus_df()
resultsDF['V (kV)'] = resultsDF['Vm'] * np.array([bus.Vnom for bus in grid9.buses])

print(resultsDF)

Which yields:

	Vm (p.u.)	V(kV)	Va (deg)	P (MW)	Q (MVAr)
Bus 1	1	17.16	0	71.63834	27.114476
Bus 2	1	18.45	9.277313	163.000003	6.973044
Bus 3	1	14.145	4.657059	85.000001	-10.669872
Bus 4	1.025709	235.91303	147.776277	-0.000009	0.000006
Bus 5	0.995498	228.96449	146.002536	-125.000009	-49.999972
Bus 6	1.01254	232.884245	146.305199	-90.000008	-29.999988
Bus 7	1.025576	235.882375	153.715955	0.000034	0.000033
Bus 8	1.015629	233.594652	150.721498	-100.000018	-34.999992
Bus 9	1.032244	237.416038	151.959026	0.000016	0.000008

Now it's on point. The yd5 vector group corresponds to 150º of phase shift which needs to be entered in rad.

Probably it'll be a good idea to implement vector groups as input.

[SanPen]

A side comment is that vector groups are often omitted for power flow but not for short circuit studies.

[fhjgch]

Thanks a lot for this, Santiago. It worked for us 👍 .

I have a follow-up question: We noticed that you chose the HELM solver. When selecting Newton-Raphson with your implementation, it suddenly doesn't converge for us. Do you know why that would be? And how could one make the example work with Newton-Raphson?

BR

[SanPen]

Hi!

As I discussed, the phase angles are usually not included in power flows because of the issue you're seeing; The solution angles became too far from the initially provided solution (Vm=1, Va=0) and iterative methods like NR cannot find it, but HELM can because it is not iterative. A kind of hacky fix is to initialize the nodal angles guess to something you estimate (see below). However it might be nicer to have the program figure that one out...

I'll see what I can do. In the mean time, this will work for you:

angles_0 = [0, 0, 0, 150, 150, 150, 150, 150, 150]
for i, bus in enumerate(grid9.buses):
    bus.Va0 = np.deg2rad(angles_0[i])

# Power flow
options = gce.PowerFlowOptions(gce.SolverType.NR,
                               retry_with_other_methods=False,
                               tolerance=1e-6,
                               control_q=False,
                               control_taps_phase=False,
                               control_taps_modules=False,
                               apply_temperature_correction=False,
                               use_stored_guess=True,
                               verbose=False)

[SanPen]

If you update to the latest version of GridCal, you'll find a new option initialize_angles which will use a linear power flow to initialize the angles, and pass on the voltages to the power flow you selected.

options = gce.PowerFlowOptions(gce.SolverType.NR,
                               retry_with_other_methods=False,
                               tolerance=1e-6,
                               control_q=False,
                               control_taps_phase=False,
                               control_taps_modules=False,
                               apply_temperature_correction=False,
                               use_stored_guess=False,
                               initialize_angles=True,
                               verbose=False)

[fhjgch]

Understood!
Thanks for the quick replies.