source: pyomo/trunk/pyomo/core/plugins/transform/radix_linearization.py @ 9475

Last change on this file since 9475 was 9475, checked in by wehart, 4 years ago

Updating the version #.

Renaming transformations to 'base.*' if they didn't have a specific
naming. NOTE: these are in the pyomo.core package ... which seems
odd. But we agreed that 'core' provides a connotation that we don't
intend.

Adding some documentation changes for 'pyomo help -t'.

File size: 11.2 KB
Line 
1#  _________________________________________________________________________
2#
3#  Pyomo: Python Optimization Modeling Objects
4#  Copyright (c) 2008 Sandia Corporation.
5#  This software is distributed under the BSD License.
6#  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
7#  the U.S. Government retains certain rights in this software.
8#  _________________________________________________________________________
9
10from pyomo.util.plugin import alias
11from pyomo.core import Binary, value, as_numeric
12from pyomo.core.base import Transformation, Var, Constraint, ConstraintList, Block, RangeSet
13from pyomo.core.base.expr import _ProductExpression, _PowExpression
14from pyomo.core.base.var import _VarData
15
16from six import iteritems
17
18import logging
19logger = logging.getLogger(__name__)
20
21class RadixLinearization(Transformation):
22    """
23    This plugin generates linear relaxations of bilinear problems using
24    the multiparametric disaggregation technique of Kolodziej, Castro,
25    and Grossmann.  See:
26
27    Scott Kolodziej, Pedro M. Castro, and Ignacio E. Grossmann. "Global
28       optimization of bilinear programs with a multiparametric
29       disaggregation technique."  J.Glob.Optim 57 pp.1039-1063. 2013.
30       (DOI 10.1007/s10898-012-0022-1)
31    """
32
33    alias( "base.radix_linearization", 
34           doc="Linearize bilinear and quadratic terms through "
35           "radix discretization (multiparametric disaggregation)" )
36
37    def apply(self, model, **kwds):
38        precision = kwds.pop('precision',8)
39        user_discretize = kwds.pop('discretize',set())
40        verbose = kwds.pop('verbose',False)
41
42        M = model.clone()
43
44        # TODO: if discretize is not empty, we must translate those
45        # components over to the components on the cloned instance
46        _discretize = {}
47        if user_discretize:
48            for _var in user_discretize:
49                _v = M.find_component(_var.cname(True))
50                if _v.component() is _v:
51                    for _vv in _v.itervalues():
52                        _discretize.setdefault(id(_vv), len(_discretize))
53                else:
54                    _discretize.setdefault(id(_v), len(_discretize))
55
56        # Iterate over all Constraints and identify the bilinear and
57        # quadratic terms
58        bilinear_terms = []
59        quadratic_terms = []
60        for constraint in M.active_components(Constraint).itervalues():
61            for cname, c in constraint._data.iteritems():
62                if c.body.polynomial_degree() != 2:
63                    continue
64                self._collect_bilinear(c.body, bilinear_terms, quadratic_terms)
65
66        # We want to find the (minimum?) number of variables to
67        # discretize so that we cover all the bilinearities -- without
68        # discretizing both sides of any single bilinear expression.
69        # First step: figure out how many expressions each term appears
70        # in
71        _counts = {}
72        for q in quadratic_terms:
73            if not q[1].is_continuous():
74                continue
75            _id = id(q[1])
76            if _id not in _counts:
77                _counts[_id] = (q[1], set())
78            _counts[_id][1].add(_id)           
79        for bi in bilinear_terms:
80            for i in (0,1):
81                if not bi[i+1].is_continuous():
82                    continue
83                _id = id(bi[i+1])
84                if _id not in _counts:
85                    _counts[_id] = (bi[i+1], set())
86                _counts[_id][1].add(id(bi[2-i]))
87
88        _tmp_counts = dict(_counts)
89        # First, remove the variables that the user wants to have discretized
90        for _id in _discretize:
91            for _i in _tmp_counts[_id][1]:
92                if _i == _id:
93                    continue
94                _tmp_counts[_i][1].remove(_id)
95            del _tmp_counts[_id]
96        # All quadratic terms must be discretized (?)
97        #for q in quadratic_terms:
98        #    _id = id(q[1])
99        #    if _id not in _tmp_counts:
100        #        continue
101        #    _discretize.setdefault(_id, len(_discretize))
102        #    for _i in _tmp_counts[_id][1]:
103        #        if _i == _id:
104        #            continue
105        #        _tmp_counts[_i][1].remove(_id)
106        #    del _tmp_counts[_id]
107           
108        # Now pick a (minimal) subset of the terms in bilinear expressions
109        while _tmp_counts:
110            _ct, _id = max( (len(_tmp_counts[i][1]), i) for i in _tmp_counts )
111            if not _ct:
112                break
113            _discretize.setdefault(_id, len(_discretize))
114            for _i in list(_tmp_counts[_id][1]):
115                if _i == _id:
116                    continue
117                _tmp_counts[_i][1].remove(_id)
118            del _tmp_counts[_id]
119
120        #
121        # Discretize things
122        #
123
124        # Define a block (namespace) for holding the disaggregated
125        # variables and new constraints
126        if False: # Set to true when the LP writer is fixed
127            M._radix_linearization = Block()
128            _block = M._radix_linearization
129        else:
130            _block = M
131        _block.DISCRETIZATION = RangeSet(precision)
132        _block.DISCRETIZED_VARIABLES = RangeSet(0, len(_discretize)-1)
133        _block.z = Var( _block.DISCRETIZED_VARIABLES, _block.DISCRETIZATION, 
134                         within=Binary )
135        _block.dv = Var( _block.DISCRETIZED_VARIABLES, 
136                         bounds=(0,2**-precision) )
137       
138        # Actually discretize the terms we have marked for discretization
139        for _id, _idx in iteritems(_discretize):
140            if verbose:
141                logger.info("Discretizing variable %s as %s" % 
142                            (_counts[_id][0].cname(True), _idx))
143            self._discretize_variable(_block, _counts[_id][0], _idx)
144
145        _known_bilinear = {}
146        # For each quadratic term, if it hasn't been discretized /
147        # generated, do so, and remember the resulting W term for later
148        # use...
149        #for _expr, _x1 in quadratic_terms:
150        #    self._discretize_term( _expr, _x1, _x1,
151        #                           _block, _discretize, _known_bilinear )
152        # For each bilinear term, if it hasn't been discretized /
153        # generated, do so, and remember the resulting W term for later
154        # use...
155        for _expr, _x1, _x2 in bilinear_terms:
156            self._discretize_term( _expr, _x1, _x2, 
157                                   _block, _discretize, _known_bilinear )
158
159        # Return the discretized instance!
160        return M
161
162
163    def _discretize_variable(self, b, v, idx):
164        _lb, _ub = v.bounds
165        if _lb is None or _ub is None:
166            raise RuntimeError("Couldn't discretize variable %s: missing "
167                               "finite lower/upper bounds." % (v.cname(True)))
168        _c = Constraint( 
169            expr= v == _lb + (_ub-_lb) * ( b.dv[idx] + 
170                sum(b.z[idx,k] * 2**-k for k in b.DISCRETIZATION) ) )
171        b.add_component("c_discr_v%s" % idx, _c)
172
173
174    def _discretize_term(self, _expr, _x1, _x2, _block, _discretize, _known_bilinear):
175        if id(_x1) in _discretize:
176            _v = _x1
177            _u = _x2
178        elif id(_x2) in _discretize:
179            _u = _x1
180            _v = _x2
181        else:
182            raise RuntimeError("Couldn't identify discretized variable "
183                               "for expression '%s'!" % _expr)
184        _id = (id(_v), id(_u))
185        if _id not in _known_bilinear:
186            _known_bilinear[_id] = self._discretize_bilinear(
187                _block, _v, _discretize[id(_v)], _u, len(_known_bilinear))
188        # _expr should be a "simple" product expression; substitute
189        # in the bilinear "W" term for the raw bilinear terms
190        _expr._numerator = [ _known_bilinear[_id] ]
191
192
193    def _discretize_bilinear(self, b, v, v_idx, u, u_idx):
194        _z = b.z
195        _dv = b.dv[v_idx]
196        _u = Var(b.DISCRETIZATION, within=u.domain, bounds=u.bounds)
197        logger.info("Discretizing (v=%s)*(u=%s) as u%s_v%s" % (
198                v.cname(True), u.cname(True), u_idx, v_idx ))
199        b.add_component( "u%s_v%s" % (u_idx, v_idx), _u)
200        _lb, _ub = u.bounds
201        if _lb is None or _ub is None:
202             raise RuntimeError("Couldn't relax variable %s: missing "
203                               "finite lower/upper bounds." % (u.cname(True)))
204        _c = ConstraintList(noruleinit=True)
205        b.add_component( "c_disaggregate_u%s_v%s" % (u_idx, v_idx), _c )
206        for k in b.DISCRETIZATION:
207            # _lb * z[v_idx,k] <= _u[k] <= _ub * z[v_idx,k]
208            _c.add(expr= _lb*_z[v_idx,k] <= _u[k] )
209            _c.add(expr= _u[k] <= _ub*_z[v_idx,k] )
210            # _lb * (1-z[v_idx,k]) <= u - _u[k] <= _ub * (1-z[v_idx,k])
211            _c.add(expr= _lb * (1-_z[v_idx,k]) <= u - _u[k] )
212            _c.add(expr= u - _u[k] <= _ub * (1-_z[v_idx,k]))
213
214        _v_lb, _v_ub = v.bounds
215        _bnd_rng = (_v_lb*_lb, _v_lb*_ub, _v_ub*_lb, _v_ub*_ub)
216        _w = Var(bounds=(min(_bnd_rng), max(_bnd_rng)))
217        b.add_component( "w%s_v%s" % (u_idx, v_idx), _w)
218
219        K = max(b.DISCRETIZATION)
220
221        _dw = Var(bounds=( min(0, _lb*2**-K, _ub*2**-K), 
222                           max(0, _lb*2**-K, _ub*2**-K) ))
223        b.add_component( "dw%s_v%s" % (u_idx, v_idx), _dw)
224
225        _c = Constraint(expr= _w == _v_lb*u + (_v_ub-_v_lb) * (
226                sum(2**-k * _u[k] for k in b.DISCRETIZATION) + _dw ) )
227        b.add_component( "c_bilinear_u%s_v%s" % (u_idx, v_idx), _c )
228
229        _c = ConstraintList(noruleinit=True)
230        b.add_component( "c_mccormick_u%s_v%s" % (u_idx, v_idx), _c )
231        # u_lb * dv <= dw <= u_ub * dv
232        _c.add(expr= _lb*_dv <= _dw )
233        _c.add(expr= _dw <= _ub*_dv )
234        # (u-u_ub)*2^-K + u_ub*dv <= dw <= (u-u_lb)*2^-K + u_lb*dv
235        _c.add(expr= (u - _ub)*2**-K + _ub*_dv <= _dw )
236        _c.add(expr= _dw <= (u - _lb)*2**-K + _lb*_dv )
237
238        return _w
239
240    def _collect_bilinear(self, expr, bilin, quad):
241        if not expr.is_expression():
242            return
243        if type(expr) is _ProductExpression:
244            if len(expr._numerator) != 2:
245                for e in expr._numerator:
246                    self._collect_bilinear(e, bilin, quad)
247                # No need to check denominator, as this is poly_degree==2
248                return
249            if not isinstance(expr._numerator[0], _VarData) or \
250                    not isinstance(expr._numerator[1], _VarData):
251                raise RuntimeError("Cannot yet handle complex subexpressions")
252            if expr._numerator[0] is expr._numerator[1]:
253                quad.append( (expr, expr._numerator[0]) )
254            else:
255                bilin.append( (expr, expr._numerator[0], expr._numerator[1]) )
256            return
257        if type(expr) is _PowExpression and value(expr._args[1]) == 2:
258            # Note: directly testing the value of the exponent above is
259            # safe: we have already verified that this expression is
260            # polynominal, so the exponent must be constant.
261            tmp = _ProductExpression()
262            tmp._numerator = [ expr._args[0], expr._args[0] ]
263            tmp._denominator = []
264            expr._args = (tmp, as_numeric(1))
265            #quad.append( (tmp, tmp._args[0]) )
266            self._collect_bilinear(tmp, bilin, quad)
267            return 
268        # All other expression types
269        for e in expr._args:
270            self._collect_bilinear(e, bilin, quad)
271
Note: See TracBrowser for help on using the repository browser.