1#  ___________________________________________________________________________
2#
3#  Pyomo: Python Optimization Modeling Objects
4#  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC
5#  Under the terms of Contract DE-NA0003525 with National Technology and
6#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
7#  rights in this software.
8#  This software is distributed under the 3-clause BSD License.
9#  ___________________________________________________________________________
10
11"""Transformation to reformulate nonlinear models with linearity induced from
12discrete variables.
13
14Ref: Grossmann, IE; Voudouris, VT; Ghattas, O. Mixed integer linear
15reformulations for some nonlinear discrete design optimization problems.
16
17"""
18
19from __future__ import division
20
21import logging
22import textwrap
23from math import fabs
24
25from pyomo.common.collections import ComponentMap, ComponentSet
26from pyomo.common.config import (ConfigBlock, ConfigValue, NonNegativeFloat,
27                                 add_docstring_list)
28from pyomo.common.modeling import unique_component_name
29from pyomo.contrib.preprocessing.util import SuppressConstantObjectiveWarning
30from pyomo.core import (Binary, Block, Constraint, Objective, Set,
31                        TransformationFactory, Var, summation, value)
32from pyomo.core.plugins.transform.hierarchy import IsomorphicTransformation
33from pyomo.gdp import Disjunct
34from pyomo.opt import TerminationCondition as tc
35from pyomo.opt import SolverFactory
36from pyomo.repn import generate_standard_repn
37
38logger = logging.getLogger('pyomo.contrib.preprocessing')
39
40
41@TransformationFactory.register('contrib.induced_linearity',
42          doc="Reformulate nonlinear constraints with induced linearity.")
43class InducedLinearity(IsomorphicTransformation):
44    """Reformulate nonlinear constraints with induced linearity.
45
46    Finds continuous variables :math:`v` where :math:`v = d_1 + d_2 + d_3`,
47    where :math:`d`'s are discrete variables. These continuous variables may
48    participate nonlinearly in other expressions, which may then be induced to
49    be linear.
50
51    The overall algorithm flow can be summarized as:
52
53    1. Detect effectively discrete variables and the constraints that
54       imply discreteness.
55    2. Determine the set of valid values for each effectively discrete variable
56    3. Find nonlinear expressions in which effectively discrete variables
57       participate.
58    4. Reformulate nonlinear expressions appropriately.
59
60    .. note:: Tasks 1 & 2 must incorporate scoping considerations (Disjuncts)
61
62    Keyword arguments below are specified for the ``apply_to`` and
63    ``create_using`` functions.
64
65    """
66
67    CONFIG = ConfigBlock("contrib.induced_linearity")
68    CONFIG.declare('equality_tolerance', ConfigValue(
69        default=1E-6,
70        domain=NonNegativeFloat,
71        description="Tolerance on equality constraints."
72    ))
73    CONFIG.declare('pruning_solver', ConfigValue(
74        default='glpk',
75        description="Solver to use when pruning possible values."
76    ))
77
78    __doc__ = add_docstring_list(__doc__, CONFIG)
79
80    def _apply_to(self, model, **kwds):
81        """Apply the transformation to the given model."""
82        config = self.CONFIG(kwds.pop('options', {}))
83        config.set_value(kwds)
84        _process_container(model, config)
85        _process_subcontainers(model, config)
86
87
88def _process_subcontainers(blk, config):
89    for disj in blk.component_data_objects(
90            Disjunct, active=True, descend_into=True):
91        _process_container(disj, config)
92        _process_subcontainers(disj, config)
93
94
95def _process_container(blk, config):
96    if not hasattr(blk, '_induced_linearity_info'):
97        blk._induced_linearity_info = Block()
98    else:
99        assert blk._induced_linearity_info.ctype == Block
100    eff_discr_vars = detect_effectively_discrete_vars(
101        blk, config.equality_tolerance)
102    # TODO will need to go through this for each disjunct, since it does
103    # not (should not) descend into Disjuncts.
104
105    # Determine the valid values for the effectively discrete variables
106    possible_var_values = determine_valid_values(blk, eff_discr_vars, config)
107
108    # Collect find bilinear expressions that can be reformulated using
109    # knowledge of effectively discrete variables
110    bilinear_map = _bilinear_expressions(blk)
111
112    # Relevant constraints are those with bilinear terms that involve
113    # effectively_discrete_vars
114    processed_pairs = ComponentSet()
115    for v1, var_values in possible_var_values.items():
116        v1_pairs = bilinear_map.get(v1, ())
117        for v2, bilinear_constrs in v1_pairs.items():
118            if (v1, v2) in processed_pairs:
119                continue
120            _process_bilinear_constraints(
121                blk, v1, v2, var_values, bilinear_constrs)
122            processed_pairs.add((v2, v1))
123            # processed_pairs.add((v1, v2))  # TODO is this necessary?
124
125
126def determine_valid_values(block, discr_var_to_constrs_map, config):
127    """Calculate valid values for each effectively discrete variable.
128
129    We need the set of possible values for the effectively discrete variable in
130    order to do the reformulations.
131
132    Right now, we select a naive approach where we look for variables in the
133    discreteness-inducing constraints. We then adjust their values and see if
134    things are stil feasible. Based on their coefficient values, we can infer a
135    set of allowable values for the effectively discrete variable.
136
137    Args:
138        block: The model or a disjunct on the model.
139
140    """
141    possible_values = ComponentMap()
142
143    for eff_discr_var, constrs in discr_var_to_constrs_map.items():
144        # get the superset of possible values by looking through the
145        # constraints
146        for constr in constrs:
147            repn = generate_standard_repn(constr.body)
148            var_coef = sum(coef for i, coef in enumerate(repn.linear_coefs)
149                           if repn.linear_vars[i] is eff_discr_var)
150            const = -(repn.constant - constr.upper) / var_coef
151            possible_vals = set((const,))
152            for i, var in enumerate(repn.linear_vars):
153                if var is eff_discr_var:
154                    continue
155                coef = -repn.linear_coefs[i] / var_coef
156                if var.is_binary():
157                    var_values = (0, coef)
158                elif var.is_integer():
159                    var_values = [v * coef for v in range(var.lb, var.ub + 1)]
160                else:
161                    raise ValueError(
162                        '%s has unacceptable variable domain: %s' %
163                        (var.name, var.domain))
164                possible_vals = set(
165                    (v1 + v2 for v1 in possible_vals for v2 in var_values))
166            old_possible_vals = possible_values.get(eff_discr_var, None)
167            if old_possible_vals is not None:
168                possible_values[eff_discr_var] = old_possible_vals & possible_vals
169            else:
170                possible_values[eff_discr_var] = possible_vals
171
172    possible_values = prune_possible_values(block, possible_values, config)
173
174    return possible_values
175
176
177def prune_possible_values(block_scope, possible_values, config):
178    # Prune the set of possible values by solving a series of feasibility
179    # problems
180    top_level_scope = block_scope.model()
181    tmp_name = unique_component_name(
182        top_level_scope, '_induced_linearity_prune_data')
183    tmp_orig_blk = Block()
184    setattr(top_level_scope, tmp_name, tmp_orig_blk)
185    tmp_orig_blk._possible_values = possible_values
186    tmp_orig_blk._possible_value_vars = list(v for v in possible_values)
187    tmp_orig_blk._tmp_block_scope = (block_scope,)
188    model = top_level_scope.clone()
189    tmp_clone_blk = getattr(model, tmp_name)
190    for obj in model.component_data_objects(Objective, active=True):
191        obj.deactivate()
192    for constr in model.component_data_objects(
193            Constraint, active=True, descend_into=(Block, Disjunct)):
194        if constr.body.polynomial_degree() not in (1, 0):
195            constr.deactivate()
196    if block_scope.ctype == Disjunct:
197        disj = tmp_clone_blk._tmp_block_scope[0]
198        disj.indicator_var.fix(1)
199        TransformationFactory('gdp.bigm').apply_to(model)
200    tmp_clone_blk.test_feasible = Constraint()
201    tmp_clone_blk._obj = Objective(expr=1)
202    for eff_discr_var, vals in tmp_clone_blk._possible_values.items():
203        val_feasible = {}
204        for val in vals:
205            tmp_clone_blk.test_feasible.set_value(eff_discr_var == val)
206            with SuppressConstantObjectiveWarning():
207                res = SolverFactory(config.pruning_solver).solve(model)
208            if res.solver.termination_condition is tc.infeasible:
209                val_feasible[val] = False
210        tmp_clone_blk._possible_values[eff_discr_var] = set(
211            v for v in tmp_clone_blk._possible_values[eff_discr_var]
212            if val_feasible.get(v, True))
213    for i, var in enumerate(tmp_orig_blk._possible_value_vars):
214        possible_values[var] = tmp_clone_blk._possible_values[
215            tmp_clone_blk._possible_value_vars[i]]
216
217    return possible_values
218
219
220def _process_bilinear_constraints(block, v1, v2, var_values, bilinear_constrs):
221    # TODO check that the appropriate variable bounds exist.
222    if not (v2.has_lb() and v2.has_ub()):
223        logger.warning(textwrap.dedent("""\
224            Attempting to transform bilinear term {v1} * {v2} using effectively
225            discrete variable {v1}, but {v2} is missing a lower or upper bound:
226            ({v2lb}, {v2ub}).
227            """.format(v1=v1, v2=v2, v2lb=v2.lb, v2ub=v2.ub)).strip())
228        return False
229    blk = Block()
230    unique_name = unique_component_name(
231        block, ("%s_%s_bilinear" % (v1.local_name, v2.local_name))
232        .replace('[', '').replace(']', ''))
233    block._induced_linearity_info.add_component(unique_name, blk)
234    # TODO think about not using floats as indices in a set
235    blk.valid_values = Set(initialize=sorted(var_values))
236    blk.x_active = Var(blk.valid_values, domain=Binary, initialize=1)
237    blk.v_increment = Var(
238        blk.valid_values, domain=v2.domain,
239        bounds=(v2.lb, v2.ub), initialize=v2.value)
240    blk.v_defn = Constraint(expr=v2 == summation(blk.v_increment))
241
242    @blk.Constraint(blk.valid_values)
243    def v_lb(blk, val):
244        return v2.lb * blk.x_active[val] <= blk.v_increment[val]
245
246    @blk.Constraint(blk.valid_values)
247    def v_ub(blk, val):
248        return blk.v_increment[val] <= v2.ub * blk.x_active[val]
249    blk.select_one_value = Constraint(expr=summation(blk.x_active) == 1)
250    # Categorize as case 1 or case 2
251    for bilinear_constr in bilinear_constrs:
252        # repn = generate_standard_repn(bilinear_constr.body)
253
254        # Case 1: no other variables besides bilinear term in constraint. v1
255        # (effectively discrete variable) is positive.
256        # if (len(repn.quadratic_vars) == 1 and len(repn.linear_vars) == 0
257        #         and repn.nonlinear_expr is None):
258        #     _reformulate_case_1(v1, v2, discrete_constr, bilinear_constr)
259
260        # NOTE: Case 1 is left unimplemented for now, because it involves some
261        # messier logic with respect to how the transformation needs to happen.
262
263        # Case 2: this is everything else, but do we want to have a special
264        # case if there are nonlinear expressions involved with the constraint?
265        pass
266        _reformulate_case_2(blk, v1, v2, bilinear_constr)
267    pass
268
269
270def _reformulate_case_2(blk, v1, v2, bilinear_constr):
271    repn = generate_standard_repn(bilinear_constr.body)
272    replace_index = next(
273        i for i, var_tup in enumerate(repn.quadratic_vars)
274        if (var_tup[0] is v1 and var_tup[1] is v2) or
275           (var_tup[0] is v2 and var_tup[1] is v1))
276    bilinear_constr.set_value((
277        bilinear_constr.lower,
278        sum(coef * repn.linear_vars[i]
279            for i, coef in enumerate(repn.linear_coefs)) +
280        repn.quadratic_coefs[replace_index] * sum(
281            val * blk.v_increment[val] for val in blk.valid_values) +
282        sum(repn.quadratic_coefs[i] * var_tup[0] * var_tup[1]
283            for i, var_tup in enumerate(repn.quadratic_vars)
284            if not i == replace_index) +
285        repn.constant +
286        zero_if_None(repn.nonlinear_expr),
287        bilinear_constr.upper
288    ))
289
290
291def zero_if_None(val):
292    return 0 if val is None else val
293
294
295def _bilinear_expressions(model):
296    # TODO for now, we look for only expressions where the bilinearities are
297    # exposed on the root level SumExpression, and thus accessible via
298    # generate_standard_repn. This will not detect exp(x*y). We require a
299    # factorization transformation to be applied beforehand in order to pick
300    # these constraints up.
301    pass
302    # Bilinear map will be stored in the format:
303    # x --> (y --> [constr1, constr2, ...], z --> [constr2, constr3])
304    bilinear_map = ComponentMap()
305    for constr in model.component_data_objects(
306            Constraint, active=True, descend_into=(Block, Disjunct)):
307        if constr.body.polynomial_degree() in (1, 0):
308            continue  # Skip trivial and linear constraints
309        repn = generate_standard_repn(constr.body)
310        for pair in repn.quadratic_vars:
311            v1, v2 = pair
312            v1_pairs = bilinear_map.get(v1, ComponentMap())
313            if v2 in v1_pairs:
314                # bilinear term has been found before. Simply add constraint to
315                # the set associated with the bilinear term.
316                v1_pairs[v2].add(constr)
317            else:
318                # We encounter the bilinear term for the first time.
319                bilinear_map[v1] = v1_pairs
320                bilinear_map[v2] = bilinear_map.get(v2, ComponentMap())
321                constraints_with_bilinear_pair = ComponentSet([constr])
322                bilinear_map[v1][v2] = constraints_with_bilinear_pair
323                bilinear_map[v2][v1] = constraints_with_bilinear_pair
324    return bilinear_map
325
326
327def detect_effectively_discrete_vars(block, equality_tolerance):
328    """Detect effectively discrete variables.
329
330    These continuous variables are the sum of discrete variables.
331
332    """
333    # Map of effectively_discrete var --> inducing constraints
334    effectively_discrete = ComponentMap()
335
336    for constr in block.component_data_objects(Constraint, active=True):
337        if constr.lower is None or constr.upper is None:
338            continue  # skip inequality constraints
339        if fabs(value(constr.lower) - value(constr.upper)
340                ) > equality_tolerance:
341            continue  # not equality constriant. Skip.
342        if constr.body.polynomial_degree() not in (1, 0):
343            continue  # skip nonlinear expressions
344        repn = generate_standard_repn(constr.body)
345        if len(repn.linear_vars) < 2:
346            # TODO should this be < 2 or < 1?
347            # TODO we should make sure that trivial equality relations are
348            # preprocessed before this, or we will end up reformulating
349            # expressions that we do not need to here.
350            continue
351        non_discrete_vars = list(v for v in repn.linear_vars
352                                 if v.is_continuous())
353        if len(non_discrete_vars) == 1:
354            # We know that this is an effectively discrete continuous
355            # variable. Add it to our identified variable list.
356            var = non_discrete_vars[0]
357            inducing_constraints = effectively_discrete.get(var, [])
358            inducing_constraints.append(constr)
359            effectively_discrete[var] = inducing_constraints
360        # TODO we should eventually also look at cases where all other
361        # non_discrete_vars are effectively_discrete_vars
362
363    return effectively_discrete
364