Automate creating Piecewise objects with a for loop

[Aziz A]

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!

[Christian Rueda Mayorga]

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)