Marginal Cost Quadratic as x2 not ^2

[Tom Adams]

Thanks for all the great work on this project. One issue I can't resolve myself unfortunately is getting quadratic expressions to work.
If I assign a quadratic marginal cost to a generator, the marginal price at the bus is simply the load x 2, not squared. E.g. in the following toy example, the marginal_cost_quadratic is $1, so I would expect a marginal price of $2500 for a load of 50MW. Instead it return $100:

import pypsa

network = pypsa.Network()

network.add("Bus", "Bus 0")

 

#add a generator at bus 0

network.add("Generator", "Gen 0",

            bus="Bus 0",

            p_nom=100,

            marginal_cost=0,

            marginal_cost_quadratic=1,

)

#add a load at bus 0 of 50MW

network.add("Load", "Load",

            bus="Bus 0",

            p_set=50,

            )

network.optimize(solver_name = 'gurobi')

print(network.buses_t.marginal_price) This return $100 ($50 x 2) as the marginal price not $2500 as I would expect. Am I missing something?

[Tom Adams]

I'll reply to myself - this is because it is the marginal cost, not the actual cost - so the marginal cost is the gradient of cost increase. That will teach me to use too simplifed an example.
However, this is a simplification of what I am really trying to do, which is to use a StorageUnit as a prioxy for network losses that go as the square of the flow across an interconnector. In this case, the amount stored really does return the flow x 2, not squared

import pypsa

network = pypsa.Network()

network.add("Bus", "Bus 0")

network.add("Bus", "Bus 1")

network.add("Line", "Line 0",

            bus0="Bus 0",

            bus1="Bus 1",

            x=0.1,

            r=0.01,

            s_nom=75)

 

#add a generator at bus 0

network.add("Generator", "Gen 0",

            bus="Bus 0",

            p_nom=100,

            marginal_cost=0,

            marginal_cost_quadratic=1,

)

#add a load at bus 1 of 50MW

network.add("Load", "Load",

            bus="Bus 1",

            p_set=50,

            )

#add a store at bus 1 to represent losses of the flow **2

network.add("StorageUnit", "Line 0 losses", bus="Bus 1", p_nom=99999,  efficiency_store = 0)  

m = network.optimize.create_model()

# Assign the inflow to the stoarge unit as a function of the flow squared

flow0 = m.variables['Line-s'].sel({'snapshot':'now','Line':'Line 0'})

loss_0 = m.variables['StorageUnit-p_store'].sel({'snapshot':'now','StorageUnit':'Line 0 losses'})

m.add_constraints(loss_0 >= (0.01 * flow0 * flow0), name="Line losses")

network.optimize.solve_model(solver_name = 'gurobi')

print(network.lines_t.p0['Line 0'].values)

print(f"Expecting {0.01 * network.lines_t.p0['Line 0'].values**2} for 0.01x flow^2")

print(f"Or {0.01 * 2 * network.lines_t.p0['Line 0'].values} for 0.01x 2 x flow")

print(network.storage_units_t.p_store['Line 0 losses'].values)