A contribution to Pyomo: linearising the multiplication of two variables

[Cristian Perfumo]

Hello everyone!

I'm working on a project that requires linearising several instances of the multiplication of two variables (e.g., AxB, AxC, CxD, DxE, etc). After realising that to linearise a multiplication you need quite a few lines of code, I wasn't going to copy and paste this like a madman, so I thought I'd come up with a function to do this. It turned out to be much harder than I anticipated, as it requires the use of function decorators. But finally I succeeded and though I'd share it here in case someone needs to linearise a product in the future.

The function uses the method described on page 85 here: http://www.aimms.com/aimms/download/manuals/aimms3om_integerprogrammingtricks.pdf

The way you use it is:

lineariseMult(model, 'A', 'B', maxA, maxB)

And it adds the variable model.A_x_B to the model. You can use this variable any time you need to use the multiplication between A and B.

I hope someone finds it useful!

Cristian

from pyomo.environ import *

from pyomo.dae import *

import numpy as np

def f2(model, j, xp):

    # we not need j, but it is passed as the index for the constraint

    return xp**2

def lineariseMult(model, var1, var2, maxVar1, maxVar2):

#Linearizing model.Temp['Tank',t]*model.Flow['Tank',t]

multName = var1 + '_x_' + var2

z1name = 'z1_' + multName

z2name = 'z2_' + multName

y1name = 'y1_' + multName

y2name = 'y2_' + multName

miny1 = - 0.5 * (maxVar1 + maxVar1)

maxy1 = 0.5 * (maxVar1 + maxVar1)

miny2 = miny1

maxy2 = maxy1 

setattr(model, z1name, Var(model.Time_Set, bounds=(0, maxy1**2)))

setattr(model, z2name, Var(model.Time_Set, bounds=(0, maxy2**2)))

setattr(model, y1name, Var(model.Time_Set, bounds=(miny1, maxy1)))

setattr(model, y2name, Var(model.Time_Set, bounds=(miny2, maxy2)))

setattr(model, multName, Var(model.Time_Set, bounds=(0, maxVar1 * maxVar2)))

setattr(model, 'lin'+z1name, Piecewise(model.Time_Set, getattr(model, z1name), getattr(model, y1name), pw_pts=list(np.arange(miny1-50, maxy1+50, 5)), pw_constr_type='EQ', f_rule=f2))

setattr(model, 'lin'+z2name, Piecewise(model.Time_Set, getattr(model, z2name), getattr(model, y2name), pw_pts=list(np.arange(miny2-50, maxy2+50, 5)), pw_constr_type='EQ', f_rule=f2))

def linDef_y1(model, t):

    return getattr(model, y1name)[t] == 0.5* (getattr(model, var1)[t] + getattr(model, var2)[t])

setattr(model, y1name+"_constraint", Constraint(model.Time_Set, rule=linDef_y1))

def linDef_y2(model, t):

    return getattr(model, y2name)[t] == 0.5* (getattr(model, var1)[t] - getattr(model, var2)[t])

setattr(model, y2name+"_constraint", Constraint(model.Time_Set, rule=linDef_y2))

def linDef_mult(model, t):

    return getattr(model, multName)[t] == getattr(model, z1name)[t] - getattr(model, z2name)[t]

setattr(model, multName+"_constraint", Constraint(model.Time_Set, rule=linDef_mult))