1

I have an optimization formulation where I have multiple decision variables, each decision variable has its own quadratic cost function term. I am planning to use a piecewise linear approximation to simplify the objective function through using the 'Piecewise' function in pyomo. I managed to that in a simple toy problem where I have a single decision variable, the issue arises when I am dealing with many decision variables. Having to write a new line for each decision variable with its own 'Piecewise' function is not viable, so I am trying to automate that with a for loop similar to the way you can do that with constraints.

Here's an example toy problem of what I'm trying to do:

import numpy as np

from pyomo.environ import *
from pyomo.core import *
from pyomo.opt import SolverFactory

def cost_function_0(x):
    return x ** 2 + 3 * x + 4

def cost_function_1(x):
    return x ** 2 + 6 * x - 2

xdata = np.linspace(-10, 10, 50)

ydata_0  = list(cost_function_0(xdata))
ydata_1  = list(cost_function_1(xdata))

xdata = list(xdata)

model = ConcreteModel()

model.N = range(2)

model.X = Var(model.N, bounds=(-10, 10))
model.Y = Var(model.N)

model.piecewise_0 = Piecewise(model.Y[0],model.X[0],
                      pw_pts=xdata,
                      pw_constr_type='EQ',
                      f_rule=ydata_0,
                      pw_repn='CC')

model.piecewise_1 = Piecewise(model.Y[1],model.X[1],
                      pw_pts=xdata,
                      pw_constr_type='EQ',
                      f_rule=ydata_1,
                      pw_repn='CC')

model.obj = Objective(expr=model.Y[0] + model.Y[1], sense=minimize)

opt = SolverFactory('glpk')
obj_val = opt.solve(model)

print('Decision variables: ', model.X[0].value, model.X[1].value)
print('Objective value: ', model.Y[0].value + model.Y[1].value)

So I am trying to replace the process of manually creating the Piecewise objects (model.piecewise_0, model.piecewise_1, ....) with an automated for loop but I got no luck so far.

Thanks in advance!

Aziz A
  • 11
  • 1

1 Answers1

0

I hope this answer helps in some way:

import numpy as np

from pyomo.environ import *
from pyomo.core import *
from pyomo.opt import SolverFactory

A  = {0:1, 1:1}
B  = {0:3, 1:6}
C  = {0:4, 1:-2}
lo = {0:-10, 1:-10}
up = {0:10, 1:10}
parts = {0:50, 1:50}

def cost_function(x,a,b,c):
    return a * x ** 2 + b * x + c


model = ConcreteModel()

model.N = range(2)

def bounds_rule(model, i):
    return (lo[i], up[i])
model.X = Var(model.N, bounds=bounds_rule)
model.Y = Var(model.N)

for idx,n in enumerate(model.N):
    xdata = np.linspace(lo[n],up[n],parts[n])
    ydata = list(cost_function(xdata,A[n],B[n],C[n]))
    xdata = list(xdata)
    
    piecewise = Piecewise(model.Y[n],model.X[n],
                          pw_pts=xdata,
                          pw_constr_type='EQ',
                          f_rule=ydata,
                          pw_repn='CC')
    model.add_component('piecewise_'+ str(idx), piecewise)

def objetive(model):
    return sum (model.Y[n] for n in model.N)
model.obj = Objective(rule = objetive, sense=minimize)

opt    = SolverFactory('glpk')
obj_val = opt.solve(model)

print('Decision variables: ', model.X[0].value, model.X[1].value)
print('Objective value: ', model.Y[0].value + model.Y[1].value)
Christian Rueda
  • 65
  • 1
  • 5