source: pyomo/trunk/pyomo/core/plugins/transform/linear_dual.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: 13.1 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#  For more information, see the FAST README.txt file.
9#  _________________________________________________________________________
10
11import itertools
12
13from pyutilib.misc import Bunch
14from pyomo.util.plugin import alias
15from pyomo.core.base import Transformation, Var, Constraint, VarList, ConstraintList, Objective, Set, maximize, minimize, NonNegativeReals, NonPositiveReals, Reals
16from pyomo.repn.canonical_repn import generate_canonical_repn
17from pyomo.repn.canonical_repn import LinearCanonicalRepn
18from pyomo.core.plugins.transform.util import process_canonical_repn
19from pyomo.bilevel import SubModel
20
21import logging
22logger = logging.getLogger('pyomo.core')
23
24
25#
26# This transformation creates a new block that
27# is the dual of the specified block.  If no block is
28# specified, then the entire model is dualized.
29# This returns a new Block object.
30#
31class LinearDual_PyomoTransformation(Transformation):
32
33    alias('base.linear_dual', doc="Dualize a linear model")
34
35    def __init__(self):
36        super(LinearDual_PyomoTransformation, self).__init__()
37
38    def apply(self, instance, **kwds):
39        options = kwds.pop('options', {})
40        bname = options.get('block',None)
41        #
42        # Iterate over the model collecting variable data,
43        # until the block is found.
44        #
45        block = None
46        if block is None:
47            block = instance
48        else:
49            for (name, data) in instance.active_components(Block).items():
50                if name == bname:
51                    block = instance
52        if block is None:
53            raise RuntimeError("Missing block: "+bname)
54        #
55        # Collect variables
56        #
57        var = {}
58        for (name, data) in block.active_components().items():
59            if isinstance(data,Var):
60                var[name] = data
61        #
62        # Generate the dual
63        #
64        instance_ = self._dualize(block)
65        #
66        # Execute the preprocessor
67        #
68        instance_.preprocess()
69        return instance_
70
71    def _dualize(self, block, unfixed=[], model=True):
72        """
73        Generate the dual of a block
74        """ 
75        #
76        # Start constructing the block
77        #
78        # Variables are constraints of block
79        # Constraints are unfixed variables of block and the parent model.
80        #
81        vnames = set()
82        for (name, data) in block.active_components(Constraint).items():
83            vnames.add((name, data.is_indexed()))
84        cnames = set(unfixed)
85        for (name, data) in block.active_components(Var).items():
86            cnames.add((name, data.is_indexed()))
87        #
88        dual = SubModel()
89        for v, is_indexed in vnames:
90            if is_indexed:
91                setattr(dual, v+'_Index', Set(dimen=None))
92                setattr(dual, v, Var(getattr(dual, v+'_Index')))
93            else:
94                setattr(dual, v, Var())
95        for cname, is_indexed in cnames:
96            if is_indexed:
97                setattr(dual, cname+'_Index', Set(dimen=None))
98                setattr(dual, cname, Constraint(getattr(dual, cname+'_Index'), noruleinit=True))
99                setattr(dual, cname+'_lower_', Var(getattr(dual, cname+'_Index')))
100                setattr(dual, cname+'_upper_', Var(getattr(dual, cname+'_Index')))
101            else:
102                setattr(dual, cname, Constraint(noruleinit=True))
103                setattr(dual, cname+'_lower_', Var())
104                setattr(dual, cname+'_upper_', Var())
105        dual.construct()
106        #
107        A = {}
108        b_coef = {}
109        c_rhs = {}
110        c_sense = {}
111        d_sense = None
112        #
113        # Collect objective
114        #
115        for (oname, odata) in block.active_components(Objective).items():
116            for ndx in odata:
117                if odata[ndx].sense == maximize:
118                    o_terms = generate_canonical_repn(-1*odata[ndx].expr, compute_values=False)
119                    d_sense = minimize
120                else:
121                    o_terms = generate_canonical_repn(odata[ndx].expr, compute_values=False)
122                    d_sense = maximize
123                for i in range(len(o_terms.variables)):
124                    c_rhs[ o_terms.variables[i].component().name, o_terms.variables[i].index() ] = o_terms.linear[i]
125            # Stop after the first objective
126            break
127        #
128        # Collect constraints
129        #
130        for (name, data) in block.active_components(Constraint).items():
131            for ndx in data:
132                con = data[ndx]
133                body_terms = generate_canonical_repn(con.body, compute_values=False)
134                lower_terms = generate_canonical_repn(con.lower, compute_values=False) if not con.lower is None else None
135                upper_terms = generate_canonical_repn(con.upper, compute_values=False) if not con.upper is None else None
136                #
137                if body_terms.constant is None:
138                    body_terms.constant = 0
139                if not lower_terms is None and not lower_terms.variables is None:
140                    raise(RuntimeError, "Error during dualization:  Constraint '%s' has a lower bound that is non-constant")
141                if not upper_terms is None and not upper_terms.variables is None:
142                    raise(RuntimeError, "Error during dualization:  Constraint '%s' has an upper bound that is non-constant")
143                #
144                for i in range(len(body_terms.variables)):
145                    varname = body_terms.variables[i].component().name
146                    varndx = body_terms.variables[i].index()
147                    A.setdefault(body_terms.variables[i].component().name, {}).setdefault(varndx,[]).append( Bunch(coef=body_terms.linear[i], var=name, ndx=ndx) )
148                   
149                #
150                if not con.equality:
151                    #
152                    # Inequality constraint
153                    #
154                    #if not (upper_terms is None or upper_terms.constant is None):
155                    if lower_terms is None or lower_terms.constant is None:
156                        #
157                        # body <= upper
158                        #
159                        v = getattr(dual,name)
160                        vardata = v.add(ndx)
161                        v.domain = NonPositiveReals
162                        b_coef[name,ndx] = upper_terms.constant - body_terms.constant
163                    #elif not (lower_terms is None or lower_terms.constant is None):
164                    elif upper_terms is None or upper_terms.constant is None:
165                        #
166                        # lower <= body
167                        #
168                        v = getattr(dual,name)
169                        vardata = v.add(ndx)
170                        v.domain = NonNegativeReals
171                        b_coef[name,ndx] = lower_terms.constant - body_terms.constant
172                    else:
173                        #
174                        # lower <= body <= upper
175                        #
176                        v = getattr(dual,name)
177                        #
178                        # Dual for lower bound
179                        #
180                        ndx_ = tuple(list(ndx).append('lb'))
181                        vardata = v.add(ndx_)
182                        vardata.domain = NonNegativeReals
183                        b_coef[name,ndx] = lower_terms.constant - body_terms.constant
184                        #
185                        # Dual for upper bound
186                        #
187                        ndx_ = tuple(list(ndx).append('ub'))
188                        vardata = v.add(ndx_)
189                        vardata.domain = NonPositiveReals
190                        b_coef[name,ndx] = upper_terms.constant - body_terms.constant
191                else:
192                    #
193                    # Equality constraint
194                    #
195                    v = getattr(dual,name)
196                    vardata = v.add(ndx)
197                    v.domain = Reals
198                    b_coef[name,ndx] = lower_terms.constant - body_terms.constant
199        #
200        # Collect bound constraints
201        #
202        for (name, data) in itertools.chain(block.active_components(Var).items(), block._parent().active_components(Var).items()):
203            #
204            # Skip fixed variables (in the parent)
205            #
206            if not (name, data.is_indexed()) in cnames:
207                continue
208            #
209            # Iterate over all variable indices
210            #
211            for ndx in data:
212                var = data[ndx]
213                bounds = var.bounds
214                if bounds[0] is None and bounds[1] is None:
215                    c_sense[name,ndx] = 'e'
216                elif bounds[0] is None:
217                    if bounds[1] == 0.0:
218                        c_sense[name,ndx] = 'g'
219                    else:
220                        c_sense[name,ndx] = 'e'
221                        #
222                        # Add constraint that defines the upper bound
223                        #
224                        name_ = name + "_upper_"
225                        varname = data.component().name
226                        varndx = data[ndx].index()
227                        A.setdefault(varname, {}).setdefault(varndx,[]).append( Bunch(coef=1.0, var=name_, ndx=ndx) )
228                        #
229                        v = getattr(dual,name_)
230                        vardata = v.add(ndx)
231                        v.domain = NonPositiveReals
232                        b_coef[name_,ndx] = bounds[1]
233                elif bounds[1] is None:
234                    if bounds[0] == 0.0:
235                        c_sense[name,ndx] = 'l'
236                    else:
237                        c_sense[name,ndx] = 'e'
238                        #
239                        # Add constraint that defines the lower bound
240                        #
241                        name_ = name + "_lower_"
242                        varname = data.component().name
243                        #from pyomo.core.base.component import Component
244                        varndx = data[ndx].index()
245                        A.setdefault(varname, {}).setdefault(varndx,[]).append( Bunch(coef=1.0, var=name_, ndx=ndx) )
246                        #
247                        v = getattr(dual,name_)
248                        vardata = v.add(ndx)
249                        v.domain = NonNegativeReals
250                        b_coef[name_,ndx] = bounds[0]
251                else:
252                    # Bounded above and below
253                    c_sense[name,ndx] = 'e'
254                    #
255                    # Add constraint that defines the upper bound
256                    #
257                    name_ = name + "_upper_"
258                    varname = data.component().name
259                    varndx = data[ndx].index()
260                    A.setdefault(varname, {}).setdefault(varndx,[]).append( Bunch(coef=1.0, var=name_, ndx=ndx) )
261                    #
262                    v = getattr(dual,name_)
263                    vardata = v.add(ndx)
264                    v.domain = NonPositiveReals
265                    b_coef[name_,ndx] = bounds[1]
266                    #
267                    # Add constraint that defines the lower bound
268                    #
269                    name_ = name + "_lower_"
270                    varname = data.component().name
271                    varndx = data[ndx].index()
272                    A.setdefault(varname, {}).setdefault(varndx,[]).append( Bunch(coef=1.0, var=name_, ndx=ndx) )
273                    #
274                    v = getattr(dual,name_)
275                    vardata = v.add(ndx)
276                    v.domain = NonNegativeReals
277                    b_coef[name_,ndx] = bounds[0]
278                    #raise IOError, "Variable bounded by (%s,%s)" % (str(bounds[0]), str(bounds[1]))
279        #
280        if d_sense == minimize:
281            dual.o = Objective(expr=sum(- b_coef[name,ndx]*getattr(dual,name)[ndx] for name,ndx in b_coef), sense=d_sense)
282        else:
283            dual.o = Objective(expr=sum(b_coef[name,ndx]*getattr(dual,name)[ndx] for name,ndx in b_coef), sense=d_sense)
284        #
285        for cname in A:
286            c = getattr(dual, cname)
287            c_index = getattr(dual, cname+"_Index") if c.is_indexed() else None
288            for ndx,terms in A[cname].iteritems():
289                if not c_index is None and not ndx in c_index:
290                    c_index.add(ndx)
291                expr = 0
292                for term in terms:
293                    v = getattr(dual,term.var)
294                    if not term.ndx in v:
295                        v.add(term.ndx)
296                    expr += term.coef * v[term.ndx]
297                if not (cname, ndx) in c_rhs:
298                    c_rhs[cname, ndx] = 0.0
299                if c_sense[cname,ndx] == 'e':
300                    c.add(ndx, expr - c_rhs[cname,ndx] == 0)
301                elif c_sense[cname,ndx] == 'l':
302                    c.add(ndx, expr - c_rhs[cname,ndx] <= 0)
303                else:
304                    c.add(ndx, expr - c_rhs[cname,ndx] >= 0)
305        return dual
306
Note: See TracBrowser for help on using the repository browser.