CODEWITHSUNDEEP
pyomo
William Bradley
William Bradley

Reputation: 11

How to Formulate a Piecewise Step Function in pyomo

I have a question regarding the correct formulation of a piecewise step function in pyomo. I want to include in my model a single piecewise function of the form:

         / 1    , 0 <= X(t) <= 1
Z(X) =   \ 0    , 1 <= X(t) <= 2

Where X is being fit to data over taken over a time domain and Z acts like a binary variable. The most similar example in pyomo documentation is the step.py example using INC. However, when solving with this formulation I observe the problem of the domain variable x ‘sticking’ to the breakpoint at x=1. I assume this is because (as noted in the documentation) Z can solve to the entire vertical line if continuous or is doubly feasible at both 0 and 1 if binary. Other formulations offered via the piecewise function (i.e. dlog, dcc, log, etc.) experience similar issues (in fact, based on the output to GAMS I’m pretty sure they don’t support binary/integer variables at all).

Is there a ‘correct’ way to formulate a piecewise function in pyomo that avoids the multiple-feasibility issue at the breakpoint, thus avoiding the domain variable converging to the breakpoint? I am using BARON with solvers cplex and ipopt, however my gut tells me this formulation issue can’t be solved by simply changing solvers.

I can also send a document illustrating my observations on why the current pyomo piecewise formulations don’t support binary variables, if it would help.

Upvotes: 1

Views: 2586

Answers (2)

Tom Roth
Tom Roth

Reputation: 2074

Here's some sample code where we try to minimise the sum of the step function Z. 

model = ConcreteModel()  
model.A = Set(initialize=[1,2,3])
model.B = Set(initialize=['J', 'K'])

model.x = Var(model.A, model.B, bounds=(0, 2))
model.z = Var(model.A, model.B, domain = Binary) 

DOMAIN_PTS = [0,1,1,2]
RANGE_PTS  = [1,1,0,0]
model.z_constraint = Piecewise(
      model.A, model.B,
      model.z, model.x,
      pw_pts=DOMAIN_PTS,
      pw_repn='INC',
      pw_constr_type = 'EQ',
      f_rule = RANGE_PTS,
      unbounded_domain_var = True)

def objective_rule(model):
    return sum(model.z[a,b] for a in model.A for b in model.B)
model.objective = Objective(rule = objective_rule, sense=minimize)

If you set sense = minimize above, the program will solve and give x = 1 for each index value. If you set sense = maximize, the program will solve and give x = 0 for each index value. I'm not too sure what you mean by stickiness, but I don't think this program does it. and it implements the step function.

Upvotes: 2

Qi Chen
Qi Chen

Reputation: 1718

This assumes that your z is not also indexed by time. If so, I would need to edit this answer:

model.t = RangeSet(*time*)
model.x = Var(model.t, bounds=(0, 2))
model.z = Var(domain=Binary)
model.d = Disjunction(expr=[
    [0 <= model.x[t] for t in model.t] + [model.x[t] <= 1 for t in model.t],
    [1 <= model.x[t] for t in model.t] + [model.x[t] <= 2 for t in model.t]
])

TransformationFactory('gdp.bigm').apply_to(model)

SolverFactory('baron').solve(model)

Upvotes: 0

Related Questions