Chain of Linear Expressions: Dynamically updating subsequent linear expressions

Dear Support-Team,

I am trying to model my problem using Gurobi (Matrix) Linear Expressions and am not sure whether the functionality I am looking for does not exist, or whether my approach is the problem. Please find the problem description below, I am stuck right now and would gratefully appreciate any help I can get.

Problem Description

My decision variable x should be maximized. In every simulation / optimization step, the constraints on x will be updated. Based on x, I am using Matrix Linear Expressions to derive values that are constrained. While these constraints remain the same in every optimization step, the way the values are calculated based on x changes. In particular, the linear expressions are calculated using x and some external data (numpy arrays) - and the values for the latter are changing with each step. Since my model is quite large, I cannot afford to re-create all linear expressions in every step. Instead, I am looking for a way to create them once and then modify them using new external data (i.e. a numpy array with new values). My issue is that I have a chain of linear expressions that depend on each other. Although I only need to modify the external values in the first linear expression (e.g. by re-creating it), the changes are not reflected in the subsequent expressions that depend on it.

More specifically, I am trying to calculate the real- and reactive power flows at nodes in a network based on the quantity traded (x) in the power market between seller i and buyer j at time t.

# basic parameters
N, T = 22, 12 # N: nodes, T: timesteps
I, J = 3, 3 # I: sellers, J: buyers
max_orders = 5 # max. possible value for I and J, respectively

# decision variable x: amount from sell order i to buy order j at time t
x = m.addMVar((max_orders, max_orders, T), lb=0, ub=M, vtype=gp.GRB.CONTINUOUS, name="x")

# objective: maximize x
m.setObjective(x.sum(), gp.GRB.MAXIMIZE)

# ... constraint x on time, price, etc. (not relevant here)

# based on where the sell- and buy orders come from, calculate the energy sold/bought at each node n at time t
energy_sold_expressions = [[gp.LinExpr() for t in range(T)] for n in range(N)]
energy_bought_expressions = [[gp.LinExpr() for t in range(T)] for n in range(N)]
for n in range(N):
    for t in range(T):
        sell_orders_at_n = [i for i in range(I) if sell_order_nodes[i] == n]
        buy_orders_at_n = [j for j in range(J) if buy_order_nodes[j] == n]
        energy_sold = gp.quicksum(x[i,j,t] for i in sell_orders_at_n for j in range(J))
        energy_bought = gp.quicksum(x[i,j,t] for i in range(I) for j in buy_orders_at_n)
        energy_sold_expressions[n][t] = energy_sold
        energy_bought_expressions[n][t] = energy_bought

# use these linear expressions to calculate the net real- and reactive power flow at each node n at time t
pf: 0.92 # pf: power factor
for n in range(N):
    for t in range(T):
        energy_sold = energy_sold_expressions[n][t]
        energy_bought = energy_bought_expressions[n][t]
        
        # calculate net real- and reactive power flow for each node and timestep
        net_real_flow = energy_sold - energy_bought
        net_reac_flow = net_real_flow * pf

# ... use net real- and reactive power flow to calculate nodal voltages, etc.

m.optimize()

In every simulation step, the amount of buy and sell orders (I, J) as well as the nodes where they come from changes (sell_order_nodes, buy_order_nodes). This means that the linear expressions that depend on them (energy_sold_expressions and energy_bought_expressions) must be updated (re-created).

However, once these variables are re-created, the new values should be used in the subsequent calculations for net real- and reactive power flows (i.e. the second loop).

Ideally, I would like to do something like this:

# new timestep, get new values for sell_order_nodes and buy_order_nodes

# re-create the first set of linear expressions
energy_sold_expressions = [[gp.LinExpr() for t in range(T)] for n in range(N)]
energy_bought_expressions = [[gp.LinExpr() for t in range(T)] for n in range(N)]
for n in range(N):
    for t in range(T):
        sell_orders_at_n = [i for i in range(I) if sell_order_nodes[i] == n]
        buy_orders_at_n = [j for j in range(J) if buy_order_nodes[j] == n]
        energy_sold = gp.quicksum(x[i,j,t] for i in sell_orders_at_n for j in range(J))
        energy_bought = gp.quicksum(x[i,j,t] for i in range(I) for j in buy_orders_at_n)
        energy_sold_expressions[n][t] = energy_sold
        energy_bought_expressions[n][t] = energy_bought

# ...somehow make sure these values are used by the subsequent linear expressions... (?)

# update and run the model again
m.optimize()

Is there any chance I can use existing Gurobi functionality to update the subsequent linear expressions without re-creating them from scratch? The model should run for 1e6+ timesteps, so this is not really an option for me. If this is not possbile, is there any change in the modelling approach that I can try?

I hope my post is understandable and would greatly appreciate any help or feedback, since I am running out of options. Thank you very much!