1

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.

William Bradley
  • 11
  • 2

2 Answers2

2

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.

Tom Roth
  • 1,954
  • 17
  • 25
  • I think the issue is the fact that I have the more complicated scenario of an objective function regressing variable X to data as part of a DAE model. Instead of smooth curve between the data points, X will converge to whatever I specify as the breakpoint in the piecewise model such that, when plotted, it generates a flat line at X=breakpoint for all indexed set t. I just wanted to know if there was a formulation that avoided this issue. – William Bradley Aug 24 '18 at 15:13
  • I've searching this for three days, you saved my life – Camunatas Nov 10 '20 at 12:35
0

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)
Qi Chen
  • 1,648
  • 1
  • 10
  • 19
  • Yes, sorry, the variable Z is indexed by time. Also, the variable X will be differentiated in another equation for dX/dt. I presume this means model.t will need to be defined by a ContinuousSet(). Are there any other changes required? – William Bradley Aug 16 '18 at 12:50
  • Yes, the z and your disjunction would also be indexed by time, and there would be a DAE transformation before solving. I don't have internet access on my computer today because I am moving, but I will make the required edits to my answer and comment again when I can. – Qi Chen Aug 17 '18 at 22:28
  • Okay. I redid the formulation to index Z by model.t = ContinuouSet(). I indexed the disjunction by adapting it to a form similar to the pyomo book example in section 9.4. (The code is too long for comments.) Unfortunately, running this code yields the following error “Found component d of type indexed by a ContinuousSet. Components of this type are not currently supported by the automatic discretization transformation in pyomo.dae. Try adding the component to the model after discretizing. Alert the pyomo developers for more assistance.” Any thoughts? – William Bradley Aug 24 '18 at 14:47
  • Hi @William, did you have any luck resolving this issue? I had a similar question – Jwem93 Oct 02 '21 at 00:12