When generator sets committable=True, buses_t.marginal_price becomes NaN

[loongmxbt]

Checklist

• I am using the current master branch or the latest release. Please indicate:
  pypsa.version
  0.16.0

Describe the Bug

Please provide a description of what the bug is and add a minimal example/command for reproducing the bug.
When there is no Generator set committable=True, everything works fine.

import pypsa

network = pypsa.Network()

network.set_snapshots(range(7))

network.add("Bus","B1")
network.add("Bus","B2")
network.add("Bus","B3")

network.add("Line", "Line1", bus0="B1", bus1="B2", x=0.01, s_nom = 150)
network.add("Line", "Line2", bus0="B1", bus1="B3", x=0.02, s_nom = 300)
network.add("Line", "Line3", bus0="B2", bus1="B3", x=0.01, s_nom = 1000)

network.add("Generator","G1", bus="B1", marginal_cost=200, p_nom=400,
#             committable=True
           )

network.add("Generator","G2", bus="B2", marginal_cost=300, p_nom=1000,
#             committable=True
           )

network.add("Load","L1", bus="B3",p_set=[280,400,401,700,701,900,901])

network.lopf()
print(network.buses)
print(network.generators_t.p)
print(network.generators_t.status)
print(network.lines_t.p0)
print(network.buses_t.marginal_price)

Below is result

attribute  v_nom type    x    y carrier  v_mag_pu_set  v_mag_pu_min  \
B1           1.0       0.0  0.0      AC           1.0           0.0   
B2           1.0       0.0  0.0      AC           1.0           0.0   
B3           1.0       0.0  0.0      AC           1.0           0.0   

attribute  v_mag_pu_max control sub_network generator  
B1                  inf   Slack           0        G1  
B2                  inf      PQ           0       NaN  
B3                  inf      PQ           0       NaN  
           G1          G2
0  280.000000    0.000000
1  333.333333   66.666667
2  333.666667   67.333333
3  400.000000  300.000000
4  400.000000  301.000000
5  300.000000  600.000000
6  299.000000  602.000000
Empty DataFrame
Columns: []
Index: [0, 1, 2, 3, 4, 5, 6]
          Line1       Line2       Line3
0  1.400000e+02  140.000000  140.000000
1  1.500000e+02  183.333333  216.666667
2  1.500000e+02  183.666667  217.333333
3  1.250000e+02  275.000000  425.000000
4  1.247500e+02  275.250000  425.750000
5 -5.684342e-13  300.000000  600.000000
6 -1.000000e+00  300.000000  601.000000
      B1     B2          B3
0  200.0  200.0  200.000000
1  200.0  300.0  266.666667
2  200.0  300.0  266.666667
3  300.0  300.0  300.000000
4  300.0  300.0  300.000000
5  200.0  300.0  400.000000
6  200.0  300.0  400.000000

When I set Generator.committable=True

import pypsa

network = pypsa.Network()

network.set_snapshots(range(7))

network.add("Bus","B1")
network.add("Bus","B2")
network.add("Bus","B3")

network.add("Line", "Line1", bus0="B1", bus1="B2", x=0.01, s_nom = 150)
network.add("Line", "Line2", bus0="B1", bus1="B3", x=0.02, s_nom = 300)
network.add("Line", "Line3", bus0="B2", bus1="B3", x=0.01, s_nom = 1000)

network.add("Generator","G1", bus="B1", marginal_cost=200, p_nom=400,
            committable=True
           )

network.add("Generator","G2", bus="B2", marginal_cost=300, p_nom=1000,
            committable=True
           )

network.add("Load","L1", bus="B3",p_set=[280,400,401,700,701,900,901])

network.lopf()
print(network.buses)
print(network.generators_t.p)
print(network.generators_t.status)
print(network.lines_t.p0)
print(network.buses_t.marginal_price)

Below is result

attribute  v_nom type    x    y carrier  v_mag_pu_set  v_mag_pu_min  \
B1           1.0       0.0  0.0      AC           1.0           0.0   
B2           1.0       0.0  0.0      AC           1.0           0.0   
B3           1.0       0.0  0.0      AC           1.0           0.0   

attribute  v_mag_pu_max control sub_network generator  
B1                  inf   Slack           0        G1  
B2                  inf      PQ           0       NaN  
B3                  inf      PQ           0       NaN  
           G1          G2
0  280.000000   -0.000000
1  333.333333   66.666667
2  333.666667   67.333333
3  400.000000  300.000000
4  400.000000  301.000000
5  300.000000  600.000000
6  299.000000  602.000000
    G1   G2
0  1.0  0.0
1  1.0  1.0
2  1.0  1.0
3  1.0  1.0
4  1.0  1.0
5  1.0  1.0
6  1.0  1.0
    Line1       Line2       Line3
0  140.00  140.000000  140.000000
1  150.00  183.333333  216.666667
2  150.00  183.666667  217.333333
3  125.00  275.000000  425.000000
4  124.75  275.250000  425.750000
5    0.00  300.000000  600.000000
6   -1.00  300.000000  601.000000
   B1  B2  B3
0 NaN NaN NaN
1 NaN NaN NaN
2 NaN NaN NaN
3 NaN NaN NaN
4 NaN NaN NaN
5 NaN NaN NaN
6 NaN NaN NaN

Error Message

If applicable, paste any terminal output to help illustrating your problem.
In some cases it may also be useful to share your list of installed packages: conda list.

print(network.buses_t.marginal_price)

   B1  B2  B3
0 NaN NaN NaN
1 NaN NaN NaN
2 NaN NaN NaN
3 NaN NaN NaN
4 NaN NaN NaN
5 NaN NaN NaN
6 NaN NaN NaN

Is there anything I miss? I'm really new to PyPSA and stuck here for a while.

[fneum]

@loongmxbt thanks for your report. There had been a similar question on the mailing list a while ago https://groups.google.com/forum/#!searchin/pypsa/shadow|sort:date/pypsa/jdeQHeZljMs/5_YUlWypCQAJ

Commitable generators make the problem mixed-integer for which shadow prices are generally not well-defined. You could determine approximate shadow prices by fixing the binary unit commitment variables. The link above has an example notebook.

import pypsa

nu = pypsa.Network()

nu.set_snapshots(range(4))

nu.add("Bus","bus")

nu.add("Generator","coal",bus="bus",
       committable=True,
       p_min_pu=0.3,
       marginal_cost=20,
       min_down_time=2,
       start_up_cost=5000,
       p_nom=10000)

nu.add("Generator","gas",bus="bus",
       committable=True,
       marginal_cost=70,
       p_min_pu=0.1,
       initial_status=0,
       shut_down_cost=25,
       p_nom=4000)

nu.add("Load","load",bus="bus",p_set=[3000,800,3000,8000])

Now build a pyomo model corresponding to the pypsa.Network:

model = pypsa.opf.network_lopf_build_model(nu, nu.snapshots)

Prepare a pyomo optimizer object, here with cbc, and solve the MILP:

opt = pypsa.opf.network_lopf_prepare_solver(nu, solver_name="cbc")
opt.solve(model).write()

Since it's a MILP cbc does not report back any dual values

model.dual.pprint()

But the generator_status has already been determined

model.generator_status.pprint()

and by fixing it, we get an LP problem, which we can solve again:

model.generator_status.fix()
nu.results = opt.solve(model)
nu.results.write()

Now all dual variables are available:

model.dual.pprint()

and we can import them back into the pypsa data-structures

pypsa.opf.extract_optimisation_results(nu, nu.snapshots)

print(nu.objective)
print(nu.generators_t.status)
print(nu.generators_t.p)
print(nu.buses_t.marginal_price)

[fneum]

Please close if this answered your questions!

[loongmxbt]

Thank you very much! Following your code solves the problem, but I still need some time to fully understand. I'll make another comment when I go through this fix.

[aleave]

I am following the suggested solution but still struggling to make sense of the result. In the toy example below, I was expecting to see the model to dispatch only gas for the first three periods and then dispatch gas + fuel oil in the last period; this because coal -which I set with a generation cost in between- has a constraint on the minimum up time of four periods and so the model should not let the solution run it for one hour only, Instead, I see the model running coal and simply ignoring the unit commitment constraint. Am I doing anything wrong?

import pypsa

nu = pypsa.Network()

nu.set_snapshots(range(4))

nu.add("Bus","bus")


nu.add("Generator","coal",bus="bus",
   committable=True,
   p_min_pu=0.1,
   marginal_cost=30,
   up_time_before=0,
   min_up_time=4,
   p_nom=1000)

nu.add("Generator","gas",bus="bus",
   marginal_cost=20,
   p_nom=1000)

nu.add("Generator","fuel_oil",bus="bus",
   marginal_cost=40,
   p_nom=1000)


nu.add("Load","load",bus="bus",p_set=[900,900,900,1800])


model = pypsa.opf.network_lopf_build_model(nu, nu.snapshots)

opt = pypsa.opf.network_lopf_prepare_solver(nu, solver_name="glpk")
opt.solve(model).write()

model.dual.pprint()
model.generator_status.pprint()

model.generator_status.fix()
nu.results = opt.solve(model)
nu.results.write()

pypsa.opf.extract_optimisation_results(nu, nu.snapshots)

results for dispatch:
print(nu.generators_t.p)
coal gas fuel_oil
0 -0.0 900.0 0.0
1 -0.0 900.0 0.0
2 -0.0 900.0 0.0
3 800.0 1000.0 0.0

and consequently marginal prices:
print(nu.buses_t.marginal_price)
bus
0 20.0
1 20.0
2 20.0
3 30.0

[fneum]

@aleave at the end of snapshots the minimum up-time is only enforced for the remaining snapshots, if the number of remaining snapshots is less than min_up_time. So if you only have four snapshots and a minimum up time of 4 hours/snapshots, no minimum up time constraints are enforced. This choice is more suited for rolling horizon optimisation, I think.

Read https://pypsa.readthedocs.io/en/latest/optimal_power_flow.html#generator-unit-commitment-constraints

[aleave]

@fneum thanks for your reply. I am still a bit confused on how to solve this in practice, even if I run the simulation for 5 snapshots with a min_up_time of 4 I see the same issue I reported. What could I try differently?