Defining a Pyomo Constraint using values from external functions/calculations

I want to solve an optimization problem using Pyomo in which the charging profile of electric vehicles is determined in a price-optimal manner while satisfying various constraints. One of these constraints should be that the line loading in the power grid where they are connected does not increase above a limit. For the modeling and simulation of the power grid I use pandapower.

So I define the time steps for which I want to create the optimized charging profiles as a pyo.Set, the price profile as a pyo.Param, and the charging profiles of the two electric vehicles each as a pyo.Var, that is, as the decision variables.

The objective function is mapped as pyo.Objective and is written to minimize the summed charging costs of both electric vehicles.

To avoid not charging at all, I have defined constraints for the minimal example below, which specify a certain minimum amount of energy that must be charged.

Up to this point the code works, the two electric vehicles charge exactly the certain amount of energy in a price-optimized way. The code so far is shown here as a minimal example:

import pandapower as pp
import pandas as pd
import pyomo.environ as pyo


net = pp.create_empty_network()
pp.create_buses(net, nr_buses=3, vn_kv=0.4)
pp.create_ext_grid(net, bus=0)
pp.create_lines(net, from_buses=[0, 1], to_buses=[1, 2], length_km=0.05, std_type="NAYY 4x50 SE")
pp.create_loads(net, buses=[1, 2], p_mw=0.011)
pp.runpp(net)


model = pyo.ConcreteModel()
model.time_steps = pyo.Set(initialize=(0, 1, 2, 3))
model.price = pyo.Param(model.time_steps, initialize=pd.Series([0.40, 0.30, 0.35, 0.45]))
model.load1_power_kw = pyo.Var(model.time_steps, within=pyo.NonNegativeReals, bounds=[0, 11], initialize=0.0)
model.load2_power_kw = pyo.Var(model.time_steps, within=pyo.NonNegativeReals, bounds=[0, 11], initialize=0.0)


def objective(m):
    costs = 0
    for ts in m.time_steps:
        cost_load1 = m.load1_power_kw[ts] * m.price[ts]
        cost_load2 = m.load2_power_kw[ts] * m.price[ts]
        costs = costs + cost_load1 + cost_load2
    return costs


model.objective = pyo.Objective(expr=objective, sense=pyo.minimize)


def pyomo_constraint_energy_load1(m):
    energy = 0
    for ts in m.time_steps:
        energy = energy + m.load1_power_kw[ts] * 0.25
    return energy >= 4


def pyomo_constraint_energy_load2(m):
    energy = 0
    for ts in m.time_steps:
        energy = energy + m.load2_power_kw[ts] * 0.25
    return energy >= 6


model.energy_constraint_load_1 = pyo.Constraint(rule=pyomo_constraint_energy_load1)
model.energy_constraint_load_2 = pyo.Constraint(rule=pyomo_constraint_energy_load2)

I then wanted to implement the already mentioned constraint of the power grid as a pyo.constraint as well, but I run into various problems here:

1. If I access the Pyomo variables directly as in the previous constraint definitions, I get the error message that an 'implicit conversion of Pyomo numeric value to float is disabled'. The Pyomo function value() is suggested as a solution.

def pyomo_constraint_pandapower_net(m, ts):
    net.load.loc[0, "p_mw"] = m.load1_power_kw[ts] / 1000
    net.load.loc[1, "p_mw"] = m.load2_power_kw[ts] / 1000
    pp.runpp(net)
    return 0, net.res_line.loading_percent.max(), 20

model.pandapower_constraint = pyo.Constraint(model.time_steps, rule=pyomo_constraint_pandapower_net)

2. However, when I use this function (no matter in which way), the constraint is not calculated correctly. If I use the .value notation

def pyomo_constraint_pandapower_net(m, ts):
    net.load.loc[0, "p_mw"] = m.load1_power_kw[ts].value / 1000
    net.load.loc[1, "p_mw"] = m.load2_power_kw[ts].value / 1000
    print(net.load.loc[0, "p_mw"])
    print(net.load.loc[1, "p_mw"])
    pp.runpp(net)
    return 0, net.res_line.loading_percent.max(), 20

model.pandapower_constraint = pyo.Constraint(model.time_steps, rule=pyomo_constraint_pandapower_net)

3. or if I use the function pyo.value

def pyomo_constraint_pandapower_net(m, ts):
    net.load.loc[0, "p_mw"] = pyo.value(m.load1_power_kw[ts]) / 1000
    net.load.loc[1, "p_mw"] = pyo.value(m.load2_power_kw[ts]) / 1000
    print(net.load.loc[0, "p_mw"])
    print(net.load.loc[1, "p_mw"])
    pp.runpp(net)
    return 0, net.res_line.loading_percent.max(), 20

model.pandapower_constraint = pyo.Constraint(model.time_steps, rule=pyomo_constraint_pandapower_net)

In both cases, no value is assigned to the loads in pandapower - the print outputs return the value 0.0 in each case.
Also, as described in this question/answer (https://stackoverflow.com/questions/73732311/using-external-functions-within-pyomo-constraints), you should not take values of variables when defining expressions in Pyomo. I don't define an expression in my minimal example, but only a constraint, but I have already tried it with expressions and got the same error messages.

I can solve the optimization problem this way, but the constraint is calculated incorrectly, so this solution is of little help to me.

results = pyo.SolverFactory("glpk").solve(model)

model.load1_power_kw.display()
model.load2_power_kw.display()
model.pandapower_constraint.display()
model.objective.display()

print(f"\nOptimal solution: {results.solver.termination_condition == pyo.TerminationCondition.optimal}")

So to summarize my problem or question again:
How can I use values from pyo.Var in an external function or pass them to an external tool like pandapower to calculate values there with the decision variables, which I can check again in a Pyomo constraint afterwards?

Thanks already for your help and best regards!

Unfortunately, I don't think you can do what you are looking to do here. I think there is a fundamental misunderstanding of how the solving construct works.

I'm not familiar with the syntax for what you are trying to do with pandapower, but it appears you are trying to use the value of a model variable to compute a constraint. This isn't possible. The value of the model constraints are unknown when the model is constructed, which is all pyomo's job is. The problem is handed off to a solver which does its magic, but does not entertain external connections.

So, you have a couple choices.... You could capture a bunch of values within some logical range and make a look-up table, or linear, or piece-wise linear line for a constraint. Or you could switch frameworks to something that chews on arbitrarily complicated functions and give up some of the performance of a good solver.