1 | from inspect import isroutine |
---|
2 | |
---|
3 | from pyomo.util.plugin import alias |
---|
4 | |
---|
5 | from pyomo.core import * |
---|
6 | from pyomo.core.base.var import SimpleVar, _VarData |
---|
7 | from pyomo.core.base.expr import _SumExpression, _ProductExpression, \ |
---|
8 | _AbsExpression, _PowExpression |
---|
9 | from pyomo.core.base.expression import _ExpressionData |
---|
10 | from pyomo.core.base.numvalue import create_name |
---|
11 | from pyomo.core.plugins.transform.util import partial |
---|
12 | from pyomo.core.plugins.transform.hierarchy import IsomorphicTransformation |
---|
13 | from pyomo.core.plugins.transform.util import collectAbstractComponents |
---|
14 | |
---|
15 | |
---|
16 | class NonNegativeTransformation(IsomorphicTransformation): |
---|
17 | """ |
---|
18 | Creates a new, equivalent model by forcing all variables to lie in |
---|
19 | the nonnegative orthant by introducing auxiliary variables. |
---|
20 | """ |
---|
21 | |
---|
22 | alias("base.nonnegative_vars", doc="Create an equivalent model in which all variables lie in the nonnegative orthant.") |
---|
23 | |
---|
24 | def __init__(self, **kwds): |
---|
25 | kwds["name"] = kwds.pop("name", "vars") |
---|
26 | super(NonNegativeTransformation, self).__init__(**kwds) |
---|
27 | |
---|
28 | self.realSets = (Reals, PositiveReals, NonNegativeReals, NegativeReals, |
---|
29 | NonPositiveReals, PercentFraction, RealSet) |
---|
30 | # Intentionally leave out Binary, Boolean, BinarySet, and BooleanSet; |
---|
31 | # we check for those explicitly |
---|
32 | self.discreteSets = (IntegerSet, Integers, PositiveIntegers, |
---|
33 | NonPositiveIntegers, NegativeIntegers, |
---|
34 | NonNegativeIntegers) |
---|
35 | |
---|
36 | |
---|
37 | def apply(self, model, **kwds): |
---|
38 | """ |
---|
39 | Force all variables to lie in the nonnegative orthant. |
---|
40 | |
---|
41 | Required arguments: |
---|
42 | model The model to transform. |
---|
43 | |
---|
44 | Optional keyword arguments: |
---|
45 | pos_suffix The suffix applied to the 'positive' component of |
---|
46 | converted variables. Default is '_plus'. |
---|
47 | neg_suffix The suffix applied to the 'positive' component of |
---|
48 | converted variables. Default is '_minus'. |
---|
49 | """ |
---|
50 | # |
---|
51 | # Optional naming schemes |
---|
52 | # |
---|
53 | pos_suffix = kwds.pop("pos_suffix", "_plus") |
---|
54 | neg_suffix = kwds.pop("neg_suffix", "_minus") |
---|
55 | # |
---|
56 | # We first perform an abstract problem transformation. Then, if model |
---|
57 | # data is available, we instantiate the new model. If not, we construct |
---|
58 | # a mapping that can later be used to populate the new model. |
---|
59 | # |
---|
60 | nonneg = model.clone() |
---|
61 | components = collectAbstractComponents(nonneg) |
---|
62 | |
---|
63 | # Map from variable base names to a {index, rule} map |
---|
64 | constraint_rules = {} |
---|
65 | |
---|
66 | # Map from variable base names to a rule defining the domains for that |
---|
67 | # variable |
---|
68 | domain_rules = {} |
---|
69 | |
---|
70 | # Map from variable base names to its set of indices |
---|
71 | var_indices = {} |
---|
72 | |
---|
73 | # Map from fully qualified variable names to replacement expressions. |
---|
74 | # For now, it is actually a map from a variable name to a closure that |
---|
75 | # must later be evaulated with a model containing the replacement |
---|
76 | # variables. |
---|
77 | var_map = {} |
---|
78 | |
---|
79 | # |
---|
80 | # Get the constraints that enforce the bounds and domains of each |
---|
81 | # variable |
---|
82 | # |
---|
83 | for var_name in components["Var"]: |
---|
84 | var = nonneg.__getattribute__(var_name) |
---|
85 | |
---|
86 | # Individual bounds and domains |
---|
87 | orig_bounds = {} |
---|
88 | orig_domain = {} |
---|
89 | |
---|
90 | # New indices |
---|
91 | indices = set() |
---|
92 | |
---|
93 | # Map from constraint names to a constraint rule. |
---|
94 | constraints = {} |
---|
95 | |
---|
96 | # Map from variable indices to a domain |
---|
97 | domains = {} |
---|
98 | |
---|
99 | for ndx in var: |
---|
100 | # Fully qualified variable name |
---|
101 | vname = create_name(str(var_name), ndx) |
---|
102 | |
---|
103 | # We convert each index to a string to avoid difficult issues |
---|
104 | # regarding appending a suffix to tuples. |
---|
105 | # |
---|
106 | # If the index is None, this casts the index to a string, |
---|
107 | # which doesn't match up with how Pyomo treats None indices |
---|
108 | # internally. Replace with "" to be consistent. |
---|
109 | if ndx is None: |
---|
110 | v_ndx = "" |
---|
111 | else: |
---|
112 | v_ndx = str(ndx) |
---|
113 | |
---|
114 | # Get the variable bounds |
---|
115 | lb = var[ndx].lb |
---|
116 | ub = var[ndx].ub |
---|
117 | if lb is not None: |
---|
118 | lb = value(lb) |
---|
119 | if ub is not None: |
---|
120 | ub = value(ub) |
---|
121 | orig_bounds[ndx] = (lb, ub) |
---|
122 | |
---|
123 | # Get the variable domain |
---|
124 | if var[ndx].domain is not None: |
---|
125 | orig_domain[ndx] = var[ndx].domain |
---|
126 | else: |
---|
127 | orig_domain[ndx] = var.domain |
---|
128 | |
---|
129 | # Determine the replacement expression. Either a new single |
---|
130 | # variable with the same attributes, or a sum of two new |
---|
131 | # variables. |
---|
132 | # |
---|
133 | # If both the bounds and domain allow for negative values, |
---|
134 | # replace the variable with the sum of nonnegative ones. |
---|
135 | |
---|
136 | bounds_neg = (orig_bounds[ndx] == (None, None) or |
---|
137 | orig_bounds[ndx][0] is None or |
---|
138 | orig_bounds[ndx][0] < 0) |
---|
139 | domain_neg = (orig_domain[ndx] is None or |
---|
140 | orig_domain[ndx].bounds()[0] is None or |
---|
141 | orig_domain[ndx].bounds()[0] < 0) |
---|
142 | if bounds_neg and domain_neg: |
---|
143 | # Make two new variables. |
---|
144 | posVarSuffix = "%s%s" % (v_ndx, pos_suffix) |
---|
145 | negVarSuffix = "%s%s" % (v_ndx, neg_suffix) |
---|
146 | |
---|
147 | new_indices = (posVarSuffix, negVarSuffix) |
---|
148 | |
---|
149 | # Replace the original variable with a sum expression |
---|
150 | expr_dict = {posVarSuffix: 1, negVarSuffix: -1} |
---|
151 | else: |
---|
152 | # Add the new index. Lie if is 'None', since Pyomo treats |
---|
153 | # 'None' specially as a key. |
---|
154 | # |
---|
155 | # More lies: don't let a blank index exist. Replace it with |
---|
156 | # '_'. I don't actually have a justification for this other |
---|
157 | # than that allowing "" as a key will eventually almost |
---|
158 | # certainly lead to a strange bug. |
---|
159 | if v_ndx is None: |
---|
160 | t_ndx = "None" |
---|
161 | elif v_ndx == "": |
---|
162 | t_ndx = "_" |
---|
163 | else: |
---|
164 | t_ndx = v_ndx |
---|
165 | |
---|
166 | new_indices = (t_ndx,) |
---|
167 | |
---|
168 | # Replace the original variable with a sum expression |
---|
169 | expr_dict = {t_ndx: 1} |
---|
170 | |
---|
171 | # Add the new indices |
---|
172 | for x in new_indices: |
---|
173 | indices.add(x) |
---|
174 | |
---|
175 | # Replace the original variable with an expression |
---|
176 | var_map[vname] = partial(self.sumRule, |
---|
177 | var_name, |
---|
178 | expr_dict) |
---|
179 | |
---|
180 | # Enforce bounds as constraints |
---|
181 | if orig_bounds[ndx] != (None, None): |
---|
182 | cname = "%s_%s" % (vname, "bounds") |
---|
183 | tmp = orig_bounds[ndx] |
---|
184 | constraints[cname] = partial( |
---|
185 | self.boundsConstraintRule, |
---|
186 | tmp[0], |
---|
187 | tmp[1], |
---|
188 | var_name, |
---|
189 | expr_dict) |
---|
190 | |
---|
191 | # Enforce the bounds of the domain as constraints |
---|
192 | if orig_domain[ndx] != None: |
---|
193 | cname = "%s_%s" % (vname, "domain_bounds") |
---|
194 | tmp = orig_domain[ndx].bounds() |
---|
195 | constraints[cname] = partial( |
---|
196 | self.boundsConstraintRule, |
---|
197 | tmp[0], |
---|
198 | tmp[1], |
---|
199 | var_name, |
---|
200 | expr_dict) |
---|
201 | |
---|
202 | # Domain will either be NonNegativeReals, NonNegativeIntegers, |
---|
203 | # or Binary. We consider Binary because some solvers may |
---|
204 | # optimize over binary variables. |
---|
205 | if isinstance(orig_domain[ndx], RealSet): |
---|
206 | for x in new_indices: |
---|
207 | domains[x] = NonNegativeReals |
---|
208 | elif isinstance(orig_domain[ndx], IntegerSet): |
---|
209 | for x in new_indices: |
---|
210 | domains[x] = NonNegativeIntegers |
---|
211 | elif isinstance(orig_domain[ndx], BooleanSet): |
---|
212 | for x in new_indices: |
---|
213 | domains[x] = Binary |
---|
214 | else: |
---|
215 | print ("Warning: domain '%s' not recognized, " + \ |
---|
216 | "defaulting to 'Reals'") % (str(var.domain)) |
---|
217 | for x in new_indices: |
---|
218 | domains[x] = Reals |
---|
219 | |
---|
220 | constraint_rules[var_name] = constraints |
---|
221 | domain_rules[var_name] = partial(self.exprMapRule, domains) |
---|
222 | var_indices[var_name] = indices |
---|
223 | |
---|
224 | # Remove all existing variables. |
---|
225 | toRemove = [] |
---|
226 | for (attr_name, attr) in nonneg.__dict__.items(): |
---|
227 | if isinstance(attr, Var): |
---|
228 | toRemove.append(attr_name) |
---|
229 | for attr_name in toRemove: |
---|
230 | nonneg.__delattr__(attr_name) |
---|
231 | |
---|
232 | # Add the sets defining the variables, then the variables |
---|
233 | for (k, v) in var_indices.items(): |
---|
234 | sname = "%s_indices" % k |
---|
235 | nonneg.__setattr__(sname, Set(initialize=v)) |
---|
236 | nonneg.__setattr__(k, Var(nonneg.__getattribute__(sname), |
---|
237 | domain = domain_rules[k], |
---|
238 | bounds = (0, None))) |
---|
239 | |
---|
240 | # Construct the model to get the variables and their indices |
---|
241 | # recognized in the model |
---|
242 | nonneg = nonneg.create() |
---|
243 | |
---|
244 | # Safe to evaluate the modifiedVars mapping |
---|
245 | for var in var_map: |
---|
246 | var_map[var] = var_map[var](nonneg) |
---|
247 | |
---|
248 | # Map from constraint base names to maps from indices to expressions |
---|
249 | constraintExprs = {} |
---|
250 | |
---|
251 | # |
---|
252 | # Convert all modified variables in all constraints in the original |
---|
253 | # problem |
---|
254 | # |
---|
255 | for conName in components["Constraint"]: |
---|
256 | con = nonneg.__getattribute__(conName) |
---|
257 | |
---|
258 | # Map from constraint indices to a corrected expression |
---|
259 | exprMap = {} |
---|
260 | |
---|
261 | for (ndx, cdata) in con._data.items(): |
---|
262 | lower = self._walk_expr(cdata.lower, var_map) |
---|
263 | body = self._walk_expr(cdata.body, var_map) |
---|
264 | upper = self._walk_expr(cdata.upper, var_map) |
---|
265 | |
---|
266 | # Lie if ndx is None. Pyomo treats 'None' indices specially. |
---|
267 | if ndx is None: |
---|
268 | ndx = "None" |
---|
269 | |
---|
270 | # Cast indices to strings, otherwise tuples ruin everything |
---|
271 | exprMap[str(ndx)] = (lower, body, upper) |
---|
272 | |
---|
273 | # Add to list of expression maps |
---|
274 | constraintExprs[conName] = exprMap |
---|
275 | |
---|
276 | # Map from constraint base names to maps from indices to expressions |
---|
277 | objectiveExprs = {} |
---|
278 | |
---|
279 | # |
---|
280 | # Convert all modified variables in all objectives in the original |
---|
281 | # problem |
---|
282 | # |
---|
283 | for objName in components["Objective"]: |
---|
284 | obj = nonneg.__getattribute__(objName) |
---|
285 | |
---|
286 | # Map from objective indices to a corrected expression |
---|
287 | exprMap = {} |
---|
288 | |
---|
289 | for (ndx, odata) in obj._data.items(): |
---|
290 | exprMap[ndx] = self._walk_expr(odata.expr, var_map) |
---|
291 | |
---|
292 | # Add to list of expression maps |
---|
293 | objectiveExprs[objName] = exprMap |
---|
294 | |
---|
295 | |
---|
296 | # Make the modified original constraints |
---|
297 | for (conName, ruleMap) in constraintExprs.items(): |
---|
298 | # Make the set of indices |
---|
299 | sname = conName + "_indices" |
---|
300 | _set = Set(initialize=ruleMap.keys()) |
---|
301 | nonneg.__setattr__(sname, _set) |
---|
302 | _set.construct() |
---|
303 | |
---|
304 | # Define the constraint |
---|
305 | _con = Constraint( nonneg.__getattribute__(sname), |
---|
306 | rule=partial(self.exprMapRule, ruleMap) ) |
---|
307 | nonneg.__setattr__(conName, _con) |
---|
308 | _con.construct() |
---|
309 | |
---|
310 | # Make the bounds constraints |
---|
311 | for (varName, ruleMap) in constraint_rules.items(): |
---|
312 | conName = varName + "_constraints" |
---|
313 | # Make the set of indices |
---|
314 | sname = conName + "_indices" |
---|
315 | _set = Set(initialize=ruleMap.keys()) |
---|
316 | nonneg.__setattr__(sname, _set) |
---|
317 | _set.construct() |
---|
318 | |
---|
319 | # Define the constraint |
---|
320 | _con = Constraint(nonneg.__getattribute__(sname), |
---|
321 | rule=partial(self.delayedExprMapRule, ruleMap)) |
---|
322 | nonneg.__setattr__(conName, _con) |
---|
323 | _con.construct() |
---|
324 | |
---|
325 | # Make the objectives |
---|
326 | for (objName, ruleMap) in objectiveExprs.items(): |
---|
327 | # Make the set of indices |
---|
328 | sname = objName + "_indices" |
---|
329 | _set = Set(initialize=ruleMap.keys()) |
---|
330 | nonneg.__setattr__(sname, _set) |
---|
331 | _set.construct() |
---|
332 | |
---|
333 | # Define the constraint |
---|
334 | _obj = Objective(nonneg.__getattribute__(sname), |
---|
335 | rule=partial(self.exprMapRule, ruleMap)) |
---|
336 | nonneg.__setattr__(objName, _obj) |
---|
337 | _obj.construct() |
---|
338 | |
---|
339 | nonneg.preprocess() |
---|
340 | return nonneg |
---|
341 | |
---|
342 | |
---|
343 | |
---|
344 | @staticmethod |
---|
345 | def _walk_expr(expr, varMap): |
---|
346 | """ |
---|
347 | Walks an expression tree, making the replacements defined in varMap |
---|
348 | """ |
---|
349 | |
---|
350 | # Attempt to replace a simple variable |
---|
351 | if isinstance(expr, SimpleVar): |
---|
352 | if expr.name in varMap: |
---|
353 | return varMap[expr.name] |
---|
354 | else: |
---|
355 | return expr |
---|
356 | |
---|
357 | # Attempt to replace an indexed variable |
---|
358 | if isinstance(expr, _VarData): |
---|
359 | if expr.name in varMap: |
---|
360 | return varMap[expr.name] |
---|
361 | else: |
---|
362 | return expr |
---|
363 | |
---|
364 | # Iterate through the numerator and denominator of a product term |
---|
365 | if isinstance(expr, _ProductExpression): |
---|
366 | i = 0 |
---|
367 | while i < len(expr._numerator): |
---|
368 | expr._numerator[i] = NonNegativeTransformation._walk_expr( |
---|
369 | expr._numerator[i], |
---|
370 | varMap) |
---|
371 | i += 1 |
---|
372 | |
---|
373 | i = 0 |
---|
374 | while i < len(expr._denominator): |
---|
375 | expr._denominator[i] = NonNegativeTransformation._walk_expr( |
---|
376 | expr.denominator[i], |
---|
377 | varMap) |
---|
378 | i += 1 |
---|
379 | |
---|
380 | # Iterate through the terms in a sum, absolute value term, or power |
---|
381 | # term |
---|
382 | if isinstance(expr, (_SumExpression, _AbsExpression, _PowExpression, |
---|
383 | _ExpressionData)): |
---|
384 | i = 0 |
---|
385 | while i < len(expr._args): |
---|
386 | expr._args[i] = NonNegativeTransformation._walk_expr( |
---|
387 | expr._args[i], |
---|
388 | varMap) |
---|
389 | i += 1 |
---|
390 | |
---|
391 | return expr |
---|
392 | |
---|
393 | @staticmethod |
---|
394 | def boundsConstraintRule(lb, ub, attr, vars, model): |
---|
395 | """ |
---|
396 | Produces 'lb < x^+ - x^- < ub' style constraints. Designed to |
---|
397 | be made a closer through functools.partial, across lb, ub, attr, |
---|
398 | and vars. vars is a {varname: coefficient} dictionary. attr is the |
---|
399 | base variable name; that is, X[1] would be referenced by |
---|
400 | |
---|
401 | model.__getattribute__('X')[1] |
---|
402 | |
---|
403 | and so attr='X', and 1 is a key of vars. |
---|
404 | |
---|
405 | """ |
---|
406 | return (lb, |
---|
407 | sum(c * model.__getattribute__(attr)[v] \ |
---|
408 | for (v,c) in vars.items()), |
---|
409 | ub) |
---|
410 | |
---|
411 | @staticmethod |
---|
412 | def noConstraint(*args): |
---|
413 | return None |
---|
414 | |
---|
415 | @staticmethod |
---|
416 | def sumRule(attr, vars, model): |
---|
417 | """ |
---|
418 | Returns a sum expression. |
---|
419 | """ |
---|
420 | return sum(c*model.__getattribute__(attr)[v] for (v, c) in vars.items()) |
---|
421 | |
---|
422 | @staticmethod |
---|
423 | def exprMapRule(ruleMap, model, ndx=None): |
---|
424 | """ Rule intended to return expressions from a lookup table """ |
---|
425 | return ruleMap[ndx] |
---|
426 | |
---|
427 | @staticmethod |
---|
428 | def delayedExprMapRule(ruleMap, model, ndx=None): |
---|
429 | """ |
---|
430 | Rule intended to return expressions from a lookup table. Each entry |
---|
431 | in the lookup table is a functor that needs to be evaluated before |
---|
432 | returning. |
---|
433 | """ |
---|
434 | return ruleMap[ndx](model) |
---|
435 | |
---|