I've tried both ways for matrix multiplication in pyomo, especially when multiplicating pyomo variables: with @ and with sum(... for index in set). It seems no difference to me regarding the result and calculation speed (for a small dataset). However, I've been told not to use @ in pyomo for some unclear reasons.
Is there any benefit to use @ over sum(... for index in set) for matrix multiplication in pyomo?
Edit: The following is my code for a reproducible example.
#Data:
import numpy as np
import pyomo.environ as pyo
gen_index = np.array([0, 0, 2, 3, 4])
bus_index = np.array([0, 1, 2, 3, 4])
nline = 6
p_init = np.array([ 20, 85, 260, 100, 300], dtype='int64')
cost_coeff = np.array([[ 0, 14, 0],
[ 0, 15, 0],
[ 0, 30, 0],
[ 0, 40, 0],
[ 0, 10, 0]], dtype='int64')
d = np.array([ 0, 300, 300, 400, 0], dtype='int64')
gen_lb = np.array([0, 0, 0, 0, 0], dtype='int64')
gen_ub = np.array([ 40, 170, 520, 200, 600], dtype='int64')
H = np.array([[1, 1, 0, 0, 0],
[0, 0, 0, 0, 0],
[0, 0, 1, 0, 0],
[0, 0, 0, 1, 0],
[0, 0, 0, 0, 1]])
M = np.array([[ 0.19391661, -0.47589472, -0.34898946, 0. , 0.15953804],
[ 0.43758813, 0.25834285, 0.18945142, 0. , 0.36001018],
[ 0.36849527, 0.21755187, 0.15953804, 0. , -0.51954822],
[ 0.19391661, 0.52410528, -0.34898946, 0. , 0.15953804],
[ 0.19391661, 0.52410528, 0.65101054, 0. , 0.15953804],
[-0.36849527, -0.21755187, -0.15953804, 0. , -0.48045178]])
f_max = np.array([400., 426., 426., 426., 426., 240.])
omega = np.array([[ 0. , -20.58723002, -0.14235327, -6.40025597,
0. ]])
#The pyomo Model:
model = pyo.ConcreteModel()
ngen = len(gen_index)
model.ng = pyo.RangeSet(0, ngen-1)
nbus = len(bus_index)
model.nb = pyo.RangeSet(0, nbus-1)
nline = 6
model.nl = pyo.RangeSet(0, nline-1)
model.p = pyo.Var(model.ng, domain=pyo.Reals, initialize=p_init)
model.omega = pyo.Var(model.nb, domain=pyo.Reals)
model.obj = pyo.Objective(expr =
cost_coeff[:,0].sum() + sum([cost_coeff[:,1][g]*model.p[g] for g in model.ng])\
+ sum([cost_coeff[:,2][g]*model.p[g] for g in model.ng])
, sense = pyo.minimize
)
model.c1b = pyo.Constraint(expr =
sum([model.p[g] for g in model.ng ]) ==
sum([d[b] for b in model.nb])
- sum([model.omega[b] for b in model.nb])
)
def c1c_lb_rule(model, g):
return gen_lb[g] <= model.p[g]
def c1c_ub_rule(model, g):
return gen_ub[g] >= model.p[g]
model.c1c_lb = pyo.Constraint(model.ng, rule=c1c_lb_rule)
model.c1c_ub = pyo.Constraint(model.ng, rule=c1c_ub_rule)
flow_expr = M@(H@model.p + model.omega - (d))
def c1d_lb_rule(model, l):
return -f_max[l] <= flow_expr[l]
def c1d_ub_rule(model, l):
return f_max[l] >= flow_expr[l]
model.c1d_lb = pyo.Constraint(model.nl, rule=c1d_lb_rule)
model.c1d_ub = pyo.Constraint(model.nl, rule=c1d_ub_rule)
#solve
solver = pyo.SolverFactory('gurobi')
obj_dict={}
infeasible=0
for i, value in (enumerate(omega)):
for j in model.nb:
model.omega[j].fix(value[j])
res = solver.solve(model)
In model.obj I'm using sum(... for index in set) instead of @ which is used in flow_exp. Both contains the variable model.p. Previously, I only used @ and it was fine as well as sum(... for index in set).
Pyomo does provide some basic support for matrix operations with NumPy classes starting in Pyomo 6.1. This is not full support for matrix operations, but rather specific support for matrix operations between Pyomo components and NumPy ndarrays (the implementation relies on the NumPy universal function API).
Currently (as of Pyomo 6.8.2), there is little difference between the matrix notation and the use of explicit loops: behind the scenes Pyomo is expanding and executing the same loops for you. I would expect that the use of the matrix operator (@) to actually be slightly slower, as Pyomo has to ensure that the Pyomo components are indexed by compatible sets, plus there will be overhead inserting the results into NumPy arrays. There are also some experimental features in development that will not be completely compatible with matrix operations with NumPy objects (the one that immediately comes to mind is the generation/use of template expressions).
As to why you may have been advised against using @ in Pyomo models, this is likely because Pyomo itself does not support a "matrix expression system" and therefore doesn't provide internal support for matrix operations. This is a feature that has been requested and would be of interest, but no one has put together a PR for it yet.