source: coopr.pyomo/stable/coopr/pyomo/base/var.py @ 3285

Last change on this file since 3285 was 3285, checked in by wehart, 10 years ago

Merged revisions 3184-3284 via svnmerge from
https://software.sandia.gov/svn/public/coopr/coopr.pyomo/trunk

........

r3188 | jwatson | 2010-10-29 08:09:35 -0600 (Fri, 29 Oct 2010) | 3 lines


Eliminating initial domain check when constructing numeric constants and the default domain (Reals) is specified.

........

r3203 | jwatson | 2010-10-29 14:46:18 -0600 (Fri, 29 Oct 2010) | 3 lines


Fixing bugs in has_discrete_variables() method of PyomoModel?.

........

r3211 | jdsiiro | 2010-11-01 14:49:11 -0600 (Mon, 01 Nov 2010) | 4 lines


bugfixes for Blocks:

  • avoid infinite loop when adding a block to a model
  • support pretty printing of user-defined components

........

r3212 | jdsiiro | 2010-11-02 15:17:13 -0600 (Tue, 02 Nov 2010) | 5 lines


  • cleaning up the management of Block._parent_block and Component.model attributes. Adding & removing blocks now updates the model attribute on all children
  • renaming Block._setattr_exec -> Block._add_component

........

r3213 | jdsiiro | 2010-11-03 10:47:06 -0600 (Wed, 03 Nov 2010) | 2 lines


Bugfix to the PyomoLogHandler? for python 2.4 compatibility

........

r3214 | jdsiiro | 2010-11-03 15:20:47 -0600 (Wed, 03 Nov 2010) | 3 lines


There is no point logging a warning when the problem is encountered
generating a logging.info message: log the warning at the info level.

........

r3215 | wehart | 2010-11-04 23:05:17 -0600 (Thu, 04 Nov 2010) | 3 lines


Adding a simple knapsack example to illustrate the difference between
a concrete and abstract model.

........

r3216 | jwatson | 2010-11-05 09:39:19 -0600 (Fri, 05 Nov 2010) | 3 lines


Fixing error diagnostic when indexing a variable with a bad index.

........

r3219 | jwatson | 2010-11-05 16:01:23 -0600 (Fri, 05 Nov 2010) | 3 lines


Supressing a validation test with NumericConstant?. If the user specifies a value, we are (now) assuming it is actually a numeric value - otherwise, the domain check significantly inflates the run-time associated with expression tree creation. This needs to be revisited in the Coopr 2.5 re-write.

........

r3226 | wehart | 2010-11-06 21:32:59 -0600 (Sat, 06 Nov 2010) | 2 lines


Setting up example, which was never converted.

........

r3233 | wehart | 2010-11-12 15:56:28 -0700 (Fri, 12 Nov 2010) | 4 lines


Migrating OS-specific functionality into coopr.os


Adding coopr.os to the dev.ini config file.

........

r3242 | wehart | 2010-11-13 01:28:57 -0700 (Sat, 13 Nov 2010) | 4 lines


Type fix.


Updating error message.

........

r3244 | wehart | 2010-11-13 10:44:26 -0700 (Sat, 13 Nov 2010) | 2 lines


Skipping OSiL writer when not defined.

........

r3246 | wehart | 2010-11-13 10:54:15 -0700 (Sat, 13 Nov 2010) | 2 lines


bug fix.

........

r3247 | wehart | 2010-11-13 11:14:01 -0700 (Sat, 13 Nov 2010) | 2 lines


Bug fix.

........

r3248 | jwatson | 2010-11-17 13:51:47 -0700 (Wed, 17 Nov 2010) | 3 lines


Interim fixes to output of quadratic terms in LP writer - more to do, but at least the basic examples now work.

........

r3254 | jwatson | 2010-11-19 13:19:19 -0700 (Fri, 19 Nov 2010) | 3 lines


Fixed bug in LP writer involving quadratic terms involving two distinct variables. Added two new quadratic examples.

........

r3257 | jwatson | 2010-11-19 13:59:35 -0700 (Fri, 19 Nov 2010) | 3 lines


Fixing diagnostic error message when attempting to solve quadratic programs with GLPK - code for generating message was not syntatically legal.

........

r3268 | jwatson | 2010-12-01 15:08:28 -0700 (Wed, 01 Dec 2010) | 3 lines


Fixing issues with the Piecewise construct when breakpoints and slopes are generated via rules. Works now (on a sample of size 1 - the newly added example5.py) for non-indexed rules, likely broken for indexed breakpoint/slope rules.

........

r3272 | jwatson | 2010-12-02 13:53:51 -0700 (Thu, 02 Dec 2010) | 3 lines


Adding omitted pprint() method for SOS constraints - identified while debugging a piecewise issue.

........

r3274 | jwatson | 2010-12-02 16:32:29 -0700 (Thu, 02 Dec 2010) | 3 lines


Adding example of Piecewise construct using breakpoint and slope rules, as opposed to explicit/direct lists.

........

r3276 | jwatson | 2010-12-03 14:06:40 -0700 (Fri, 03 Dec 2010) | 3 lines


Some progress toward functional indexed Piecewise components.

........

File size: 44.2 KB
Line 
1#  _________________________________________________________________________
2#
3#  Coopr: A COmmon Optimization Python Repository
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 Coopr README.txt file.
9#  _________________________________________________________________________
10
11__all__ = ['Var', 'VarStatus', '_VarBase', "Piecewise"]
12
13from component import Component
14from numvalue import *
15from numvalue import create_name
16from pyutilib.component.core import alias
17import pyomo
18import types
19from set_types import *
20from indexed_component import IndexedComponent
21import pyutilib.math
22import pyutilib.misc
23from pyutilib.enum import Enum
24import sys
25from param import _ParamValue
26from coopr.pyomo.base.util import isfunctor
27
28VarStatus = Enum( 'undefined', 'fixed_by_presolve', 'fixed_by_optimizer',
29                  'unused', 'used', )
30
31class _VarValue(NumericValue):
32    """Holds the numeric value of a variable, along with suffix info"""
33
34    def __init__(self,**kwds):
35        """Constructor"""
36        if kwds.pop('_deep_copying', None):
37            # Hack for Python 2.4 compatibility
38            # Deep copy will copy all items as necessary, so no need to
39            # complete parsing
40            return
41
42        # NOTE: we currently use the name (in the base class), id, and label attributes,
43        #       but don't document their differences anywhere!
44        # NOTE: the "within" keyword takes precedence over the "domain" keyword.
45        NumericValue.__init__(self, \
46                              name=kwds.get('name', None), \
47                              domain=kwds.get('within', kwds.get('domain', Reals)))
48
49        self.var = None # the "parent" variable.
50        self.index = None # the index of this variable within the "parent"
51
52        # we currently have too many names/labels floating around - the
53        # two below and the name in the base numeric value class.
54        self.id = None
55        self.label = None
56
57        self.initial = None
58
59        self.active = True
60
61        # IMPORTANT: in contrast to the initial value, the lower and upper
62        #             bounds of a variable can in principle be either numeric constants,
63        #             parameter values, or expressions - anything is OK, as long as
64        #             invoking () is legal.
65        self.lb = None
66        self.ub = None
67
68        self.fixed = False
69        self.status = VarStatus.undefined
70
71        # cached for efficiency purposes - isinstance is not cheap.
72        self._is_binary = isinstance(self.domain, BooleanSet)
73        self._is_integer = isinstance(self.domain, IntegerSet)
74        self._is_continuous = not (self._is_binary or self._is_integer)
75
76    def __str__(self):
77        return self.name
78
79    def activate(self):
80        self.active=True
81
82    def deactivate(self):
83        self.active=False
84
85    def setlb(self, value):
86        # python only has 3 numeric built-in type that we deal with (we ignore
87        # complex numbers).
88        if isinstance(value, (int, long, float)):
89            self.lb = NumericConstant(value=value)
90        elif value is None:
91            self.lb = None
92        elif isinstance(value, _ParamValue):
93            self.lb = value
94        else:
95            msg = "Unknown type '%s' supplied as variable lower bound - "     \
96                  'legal types are numeric constants or parameters'
97            raise ValueError, msg % str( type(value) )
98
99    def setub(self, value):
100        # python only has 3 numeric built-in type that we deal with (we ignore
101        # complex numbers).
102        if isinstance(value, (int, long, float)):
103            self.ub = NumericConstant(value=value)
104        elif value is None:
105            self.ub = None
106        elif isinstance(value, _ParamValue):
107            self.ub = value
108        else:
109            msg = "Unknown type '%s' supplied as variable lower bound - "     \
110                  'legal types are numeric constants or parameters'
111            raise ValueError, msg % str( type(value) )
112
113
114    def fixed_value(self):
115        if self.fixed:
116            return True
117        return False
118
119    def is_constant(self):
120        if self.fixed:
121            return True
122        return False
123
124    def is_binary(self):
125        return self._is_binary
126
127    def is_integer(self):
128        return self._is_integer
129
130    def is_continuous(self):
131        return self._is_continuous
132
133    def simplify(self, model):
134        return self                 #pragma:nocover
135
136    def pprint(self, ostream=None):
137        if ostream is None:
138           ostream = sys.stdout
139        print >>ostream, str(self),
140
141    def _tighten_bounds(self):
142        """
143        Attempts to tighten the lower- and upper-bounds on this variable
144        by using the domain.
145        """
146        try:
147            dbounds = self.domain.bounds()
148        except:
149            # No explicit domain bounds
150            return
151
152        if self.lb is None or \
153               (dbounds[0] is not None and dbounds[0] > self.lb.value):
154            self.setlb(dbounds[0])
155
156        if self.ub is None or \
157               (dbounds[1] is not None and dbounds[1] < self.ub.value):
158            self.setub(dbounds[1])
159
160
161#
162# Variable attributes:
163#
164# Single variable:
165#
166# x = Var()
167# x.fixed = True
168# x.domain = Reals
169# x.name = "x"
170# x.setlb(-1.0)
171# x.setub(1.0)
172# x.initial = 0.0
173# x = -0.4
174#
175# Array of variables:
176#
177# y = Var(set1,set2,...,setn)
178# y.fixed = True (fixes all variables)
179# y.domain = Reals
180# y.name = "y"
181# y[i,j,...,k] (returns value of a given index)
182#
183# y[i,j,...,k].fixed
184# y[i,j,...,k].lb()
185# y[i,j,...,k].ub()
186# y[i,j,...,k].initial
187#
188
189class _VarBase(IndexedComponent):
190    """A numeric variable, which may be defined over a index"""
191
192    """ Constructor
193        Arguments:
194           name         The name of this variable
195           index        The index set that defines the distinct variables.
196                          By default, this is None, indicating that there
197                          is a single variable.
198           domain       A set that defines the type of values that
199                          each parameter must be.
200           default      A set that defines default values for this
201                          variable.
202           bounds       A rule for defining bounds values for this
203                          variable.
204           rule         A rule for setting up this parameter with
205                          existing model data
206    """
207    # TODO: default and rule keywords are not used?  Need to talk to Bill ...?
208    def __init__(self, *args, **kwd):
209        if kwd.pop('_deep_copying', None):
210            # Hack for Python 2.4 compatibility
211            # Deep copy will copy all items as necessary, so no need to
212            # complete parsing
213            return
214
215        tkwd = {'ctype':Var}
216        IndexedComponent.__init__(self, *args, **tkwd)
217
218        # Default keyword values
219        #
220        self._attr_declarations = {}
221        self._varval = {}
222        self._initialize = kwd.pop('initialize', None )
223        self.name   = kwd.pop('name', 'unknown')
224        self.domain = kwd.pop('within', Reals )
225        self.domain = kwd.pop('domain', self.domain )
226        self.bounds = kwd.pop('bounds', None )
227        self.doc    = kwd.pop('doc', None )
228
229        self._binary_keys = []
230        self._continuous_keys = []
231        self._integer_keys = []
232
233        # Check for domain rules
234        if isfunctor(self.domain):
235            self._domain_rule = self.domain
236
237            # Individual variables will be restricted
238            self.domain = None
239        else:
240            self._domain_rule = None
241
242        if kwd:
243            msg = "Var constructor: unknown keyword '%s'"
244            raise ValueError, msg % kwd.keys()[0]
245
246    def as_numeric(self):
247        if None in self._varval:
248            return self._varval[None]
249        return self
250
251    def keys(self):
252        return self._varval.keys()
253
254    def binary_keys(self):
255        """ Returns the keys of all binary variables """
256        return self._binary_keys
257
258    def continuous_keys(self):
259        """ Returns the keys of all continuous variables """
260        return self._continuous_keys
261
262    def integer_keys(self):
263        """ Returns the keys of all integer variables """
264        return self._integer_keys
265
266    def __iter__(self):
267        return self._varval.keys().__iter__()
268
269    def __contains__(self,ndx):
270        return ndx in self._varval
271
272    def dim(self):
273        return self._ndim
274
275    def reset(self):
276        for key in self._varval:
277            if not self._varval[key].initial is None:
278                self._varval[key].set_value( self._varval[key].initial )
279
280    def __len__(self):
281        return len(self._varval)
282
283    def __setitem__(self,ndx,val):
284        #print "HERE",ndx,val, self._valid_value(val,False), self.domain
285        if (None in self._index) and (ndx is not None):
286            # allow for None indexing if this is truly a singleton
287            msg = "Cannot set an array value in singleton variable '%s'"
288            raise KeyError, msg % self.name
289
290        if ndx not in self._index:
291            msg = "Cannot set the value of array variable '%s' with invalid " \
292                  "index '%s'"
293            raise KeyError, msg % ( self.name, str(ndx) )
294
295        if not self._valid_value(val,False):
296            msg = "Cannot set variable '%s[%s]' with invalid value: %s"
297            raise ValueError, msg % ( self.name, str(ndx), str(val) )
298        self._varval[ndx].value = val
299
300    def __getitem__(self,ndx):
301        """This method returns a _VarValue object.  This object can be
302           coerced to a numeric value using the value() function, or using
303           explicity coercion with float().
304        """
305        try:
306            return self._varval[ndx]
307        except KeyError: # thrown if the supplied index is hashable, but not defined.
308            msg = "Unknown index '%s' in variable %s;" % (str(ndx), self.name)
309            if (isinstance(ndx, (tuple, list)) and len(ndx) != self.dim()):
310                msg += "    Expecting %i-dimensional indices" % self.dim()
311            else:
312                msg += "    Make sure the correct index sets were used.\n"
313                msg += "    Is the ordering of the indices correct?"
314            raise KeyError, msg
315        except TypeError, msg: # thrown if the supplied index is not hashable
316            msg2 = "Unable to index variable %s using supplied index with " % self.name
317            msg2 += str(msg)
318            raise TypeError, msg2
319     
320
321    # TODO: Convert to 'has_XXX'
322#     def is_binary(self):
323#         return isinstance(self.domain, BooleanSet)
324
325#     def is_integer(self):
326#         return isinstance(self.domain, IntegerSet)
327
328#     def is_continuous(self):
329#         return not (self.is_binary() or self.is_integer())
330
331    def simplify(self, model):
332        return self
333
334    def _add_domain_key(self, ndx, domain):
335        """ Register an index with a specific set of keys """
336        if isinstance(domain, BooleanSet):
337            self._binary_keys.append(ndx)
338        elif isinstance(domain, IntegerSet):
339            self._integer_keys.append(ndx)
340        else:
341            self._continuous_keys.append(ndx)
342
343    def construct(self, data=None):
344        if pyomo.debug("verbose"):      #pragma:nocover
345           print "Constructing Variable, name="+self.name+", from data="+`data`
346        if self._constructed:
347            return
348        self._constructed=True
349        #
350        # Construct _VarValue() objects for all index values
351        #
352        if self._ndim > 0:
353            if type(self._initialize) is dict:
354                self._index = self._initialize.keys()
355
356            for ndx in self._index:
357                # Check for domain rules
358                if self._domain_rule is not None:
359                    domain = self._domain_rule(ndx, self.model)
360                    if isinstance(domain, BooleanSet):
361                        self._binary_keys.append(ndx)
362                    elif isinstance(domain, IntegerSet):
363                        self._integer_keys.append(ndx)
364                    else:
365                        self._continuous_keys.append(ndx)
366                else:
367                    domain = self.domain
368                    self._add_domain_key(ndx, domain)
369
370                self._varval[ndx] = _VarValue(
371                      name = create_name(self.name,ndx),
372                      domain = domain,
373                      )
374
375                self._varval[ndx].var = self
376                self._varval[ndx].index = ndx
377
378        #
379        # Define the _XXX_keys objects if domain isn't a rule;
380        # they were defined individually above in the case of
381        # a rule
382        #
383        if not self._domain_rule is not None:
384            if isinstance(self.domain, BooleanSet):
385                self._binary_keys = self._varval.keys()
386            elif isinstance(self.domain, IntegerSet):
387                self._integer_keys = self._varval.keys()
388            else:
389                self._continuous_keys = self._varval.keys()
390
391        #
392        # Initialize values with a dictionary if provided
393        #
394        if self._initialize is not None and \
395           type(self._initialize) is not types.FunctionType:
396
397           if type(self._initialize) is dict:
398              for key in self._initialize:
399                self._varval[key].value = None
400                self._varval[key].initial = self._initialize[key]
401                self._valid_value(self._initialize[key],True)
402           else:
403              for key in self._index:
404                self._varval[key].value = None
405                self._varval[key].initial = self._initialize
406              self._valid_value(self._initialize,True)
407        #
408        # Initialize values with the _rule function if provided
409        #
410        elif type(self._initialize) is types.FunctionType:
411           for key in self._varval:
412             if isinstance(key,tuple):
413                tmp = list(key)
414             elif key is None:
415                tmp = []
416             else:
417                tmp = [key]
418             tmp.append(self.model)
419             tmp = tuple(tmp)
420             self._varval[key].value = None
421             self._varval[key].initial = self._initialize(*tmp)
422             self._valid_value(self._varval[key].initial,True)
423        #
424        # Initialize bounds with the bounds function if provided
425        #
426        if self.bounds is not None:
427           # bounds are specified via a tuple
428           if type(self.bounds) is tuple:
429             for key in self._varval:
430               (lb, ub) = self.bounds
431               self._varval[key].setlb(lb)
432               self._varval[key].setub(ub)
433           else:
434             # bounds are specified via a function
435             for key in self._varval:
436               if isinstance(key,tuple):
437                  tmp = list(key)
438               elif key is None:
439                  tmp = []
440               else:
441                  tmp = [key]
442               tmp.append(self.model)
443               tmp = tuple(tmp)
444               (lb, ub) = self.bounds(*tmp)
445               self._varval[key].setlb(lb)
446               self._varval[key].setub(ub)
447           for key in self._varval:
448             if self._varval[key].lb is not None and \
449                not pyutilib.math.is_finite(self._varval[key].lb()):
450                self._varval[key].setlb(None)
451             if self._varval[key].ub is not None and \
452                not pyutilib.math.is_finite(self._varval[key].ub()):
453                self._varval[key].setub(None)
454        #
455        # Iterate through all variables, and tighten the bounds based on
456        # the domain bounds information.
457        #
458        # Only done if self.domain is not a rule. If it is, _VarArray level
459        # bounds become meaningless, since the individual _VarElement objects
460        # likely have more restricted domains.
461        #
462        if self._domain_rule is None:
463            dbounds = self.domain.bounds()
464            if not dbounds is None and dbounds != (None,None):
465                for key in self._varval:
466                    if not dbounds[0] is None:
467                        if self._varval[key].lb is None or                     \
468                            dbounds[0] > self._varval[key].lb():
469                            self._varval[key].setlb(dbounds[0])
470
471                    if not dbounds[1] is None:
472                        if self._varval[key].ub is None or                     \
473                            dbounds[1] < self._varval[key].ub():
474                            self._varval[key].setub(dbounds[1])
475        #
476        # Setup declared attributes for all variables
477        #
478        for attr in self._attr_declarations:
479            for key in self._varval:
480                setattr(self._varval[key],attr,self._attr_declarations[attr][0])
481
482    def pprint(self, ostream=None):
483        if ostream is None:
484           ostream = sys.stdout
485        print >>ostream, "  ",self.name,":",
486        print >>ostream, "\tSize="+str(len(self)),
487        if self.domain is not None:
488            print >>ostream, "\tDomain="+self.domain.name
489        else:
490            print >>ostream, "\tDomain=None"
491        if self._index_set is not None:
492            print >>ostream, "\tIndicies: ",
493            for idx in self._index_set:
494               print >>ostream, str(idx.name)+", ",
495            print ""
496        if None in self._varval:
497           print >>ostream, "\tInitial Value : Lower Bound : Upper Bound : "  \
498                            "Current Value: Fixed: Status"
499           lb_value = None
500           if self._varval[None].lb is not None:
501              lb_value = self._varval[None].lb()
502           ub_value = None
503           if self._varval[None].ub is not None:
504              ub_value = self._varval[None].ub()
505
506           print >>ostream, "\t %s : %s : %s : %s : %s : %s" % (
507             str(self._varval[None].initial),
508             str(lb_value),
509             str(ub_value),
510             str(self._varval[None].value),
511             str(self._varval[None].fixed),
512             str(self._varval[None].status)
513           )
514        else:
515           print >>ostream, "\tKey : Initial Value : Lower Bound : "          \
516                            "Upper Bound : Current Value: Fixed: Status"
517           tmp=self._varval.keys()
518           tmp.sort()
519           for key in tmp:
520             initial_val = self._varval[key].initial
521             lb_value = None
522             if self._varval[key].lb is not None:
523                lb_value = self._varval[key].lb()
524             ub_value = None
525             if self._varval[key].ub is not None:
526                ub_value = self._varval[key].ub()
527
528             print >>ostream, "\t%s : %s : %s : %s : %s : %s : %s" % (
529               str(key),
530               str(initial_val),
531               str(lb_value),
532               str(ub_value),
533               str(value(self._varval[key].value)),
534               str(value(self._varval[key].fixed)),
535               str(self._varval[key].status)
536             )
537
538
539    def display(self, prefix="", ostream=None):
540        if ostream is None:
541           ostream = sys.stdout
542        print >>ostream, prefix+"Variable "+self.name,":",
543        print >>ostream, "  Size="+str(len(self)),
544        print >>ostream, "Domain="+self.domain.name
545        if None in self._varval:
546           print >>ostream, "%s  Value=%s" % (
547             prefix,
548             pyutilib.misc.format_io(self._varval[None].value)
549           )
550        else:
551           for key in self._varval:
552             val = self._varval[key].value
553             print >>ostream, prefix+"  "+str(key)+" : "+str(val)
554
555    def declare_attribute(self, name, default=None):
556        """
557        Declare a user-defined attribute.
558        """
559        if name[0] == "_":
560            msg = "Cannot define an attribute that begins with  '_'"
561            raise AttributeError, msg
562        if name in self._attr_declarations:
563            raise AttributeError, "Attribute %s is already defined" % name
564        self._attr_declarations[name] = (default,)
565        #if not default is None:
566            #self._valid_value(default)
567        #
568        # If this variable has been constructed, then
569        # generate this attribute for all variables
570        #
571        if len(self._varval) > 0:
572            for key in self._varval:
573                setattr(self._varval[key],name,default)
574
575
576class _VarElement(_VarBase,_VarValue):
577
578    def __init__(self, *args, **kwd):
579        if kwd.pop('_deep_copying', None):
580            # Hack for Python 2.4 compatibility
581            # Deep copy will copy all items as necessary, so no need to
582            # complete parsing
583            return
584
585        _VarValue.__init__(self, **kwd)
586        _VarBase.__init__(self, *args, **kwd)
587        self._varval[None] = self
588        self._varval[None].var = self
589        self._varval[None].index = None
590
591    def __call__(self, exception=True):
592        if None in self._varval:
593            return self._varval[None].value
594        return None
595
596    def is_constant(self):
597        return _VarValue.is_constant(self)
598
599class _VarArray(_VarBase):
600
601    def __init__(self, *args, **kwd):
602        if kwd.pop('_deep_copying', None):
603            # Hack for Python 2.4 compatibility
604            # Deep copy will copy all items as necessary, so no need to
605            # complete parsing
606            return
607
608        _VarBase.__init__(self, *args, **kwd)
609        self._dummy_val = _VarValue(**kwd)
610
611    def __float__(self):
612        raise TypeError, "Cannot access the value of array variable "+self.name
613
614    def __int__(self):
615        raise TypeError, "Cannot access the value of array variable "+self.name
616
617    def _valid_value(self,value,use_exception=True):
618        #print "VALID",self._dummy_val.domain
619        return self._dummy_val._valid_value(value,use_exception)
620
621    def set_value(self, value):
622        msg = "Cannot specify the value of array variable '%s'"
623        raise ValueError, msg % self.name
624
625    def __str__(self):
626        return self.name
627
628    def construct(self, data=None):
629        _VarBase.construct(self, data)
630
631        for (ndx, var) in self._varval.iteritems():
632            var._tighten_bounds()
633
634
635class Var(Component):
636    """
637    Variable objects that are used to construct Pyomo models.
638    """
639
640    alias("Var", "Decision variables in a model.")
641
642    def __new__(cls, *args, **kwds):
643        if kwds.pop('_deep_copying', None):
644            # Hack for Python 2.4 compatibility
645            # Deep copy will copy all items as necessary, so no need to
646            # complete parsing
647            return
648
649        if args == ():
650            self = _VarElement(*args, **kwds)
651        else:
652            self = _VarArray(*args, **kwds)
653        return self
654
655
656class Piecewise(_VarElement):
657    """
658    A piecewise linear expression.
659
660    Usage:
661
662    model.X = Var()
663
664    # Let breakpoint_rule and slope_rule define a piecewise linear
665    # function f(x) when indexed over sets (...)
666    #
667    # Then model.XF is a variable constrained by model.XF = f(model.X)
668    model.XF = Piecewise(
669                         breakpoint_rule=breakPointRule,
670                         slope_rule=slopeRule,
671                         breakpoint_sets=(...),
672                         slope_sets=(...),
673                         offset=0,
674                         value=model.X,
675                         forcesos2=True
676                         )
677
678    # Evaluate the piecewise expression at 7.0
679    model.C2.eval(7.0)
680
681    Breakpoints refer to the values on the domain that separate the
682    domain into distinct intervals. All intervals are assumed to be
683    closed on the lower bound and open on the upper bound, with the
684    exception of the interval extending to negative infinity on the
685    lower bound. The first slope specified is assumed to be the slope
686    of an interval coming from negative infinity to the first
687    breakpoint; all other slopes are assumed to correspond to the
688    intervals between the current and next breakpoint. The offset
689    keyword is used to indicate the y-intercept of the expression
690    (default is zero).
691
692    Example:
693
694    The concave function
695
696           / 2x + 3, x < 0
697    f(x) = | 3,      0 < x < 1
698           \ 3 - x,  1 < x
699
700    would be modeled as:
701
702    breakpoints = [0, 1]
703    slopes = [2, 0, -1]
704    offset = 3
705
706    If the function should have a bounded domain, specify more
707    breakpoints than slopes.
708
709    Example:
710
711    The concave function
712
713           / 2x + 3, -2 < x < 0
714    f(x) = | 3,       0 < x < 1
715           \ 3 - x,   1 < x < 4
716
717    would be modeled as:
718
719    breakpoints = [-2, 0, 1, 4]
720    slopes = [2, 0, -1]
721    offset = 3
722
723    If we model a bounded domain, the variable passed to model.value
724    will be bounded by this domain. Finally, for bounded domains, the
725    offset is the vertical offset of the leftmost breakpoint.
726
727    The forcesos2 keyword is simply used to always generate the sos2
728    transformation, even if the slopes specify a concave/convex function.
729
730    """
731
732    alias("Piecewise", "A variable constrained by a piecewise function.")
733
734    def __init__(self, *args, **kwargs):
735        """
736        Create a new Piecewise variable
737        """
738
739        # ---------------------------------------------------------------------
740        # Process arguments
741
742        # We can either get explicit breakpoints and slopes, or rules
743        # defining them.
744        self._breakpoints = kwargs.pop('breakpoints', None)
745        self._slopes = kwargs.pop('slopes', None)
746
747        self._breakpoint_rule = kwargs.pop('breakpoint_rule', None)
748        self._slope_rule = kwargs.pop('slope_rule', None)
749
750        _breakpoint_sets = kwargs.pop('breakpoint_sets', None)
751        _slope_sets = kwargs.pop('slope_sets', None)
752
753        self._forcesos2 = kwargs.pop('forcesos2', False)
754
755        # The corresponding domain variable.
756        self.variable = kwargs.pop('value', None)
757
758        # The y-intercept of the entire piecewise construct.
759        self.offset = kwargs.pop('offset', 0)
760
761        # Don't pop the name, we want to pass this to superclass constructor
762        tmpname = kwargs.get('name', 'unknown')
763
764        # ---------------------------------------------------------------------
765        # Error checking
766
767        # Make sure a domain variable was passed.
768        if self.variable is None:
769            raise TypeError, "Specify a Variable object for the 'value' " \
770                  "attribute for Piecewise object '%s'" % tmpname
771
772        # Make sure repetitive arguments aren't specified
773        if self._breakpoints is not None and self._breakpoint_rule is not None:
774            raise TypeError, "Specify only one of 'breakpoints' and " \
775                  "'breakpointrule' in Piecewise object '%s'" % tmpname
776        if self._slopes is not None and self._slope_rule is not None:
777            raise TypeError, "Specify only one of 'slopes' and " \
778                  "'sloperule' in Piecewise object '%s'" % tmpname
779
780        # Make sure enough arguments were specified
781        if self._breakpoints is None and self._breakpoint_rule is None:
782            raise TypeError, "Specify exactly one of 'breakpoints' or " \
783                  "'breakpointrule' in Piecewise object '%s'" % tmpname
784        if self._slopes is None and self._slope_rule is None:
785            raise TypeError, "Specify exactly one of 'slopes' or " \
786                  "'sloperule' in Piecewise object '%s'" % tmpname
787
788        # ---------------------------------------------------------------------
789        # Initialization
790
791        # Call superclass
792        _VarElement.__init__(self, *args, **kwargs)
793
794        # Create indexed_components for the breakpoints and slopes
795        if _breakpoint_sets is not None:
796            self._breakpoint_component = IndexedComponent(
797                _breakpoint_sets,
798                ctype="IndexedComponent"
799                )
800        else:
801            self._breakpoint_component = None
802        if _slope_sets is not None:
803            self._slope_component = IndexedComponent(
804                _slope_sets,
805                ctype="IndexedComponent"
806                )
807        else:
808            self._slope_component = None
809
810        # If index sets were passed to either breakpoint or slopes,
811        # make sure the appropriate rule is defined
812        if self._breakpoint_component is not None and \
813               self._breakpoint_rule is None:
814            raise TypeError, "Cannot specify 'breakpoint_sets' without " \
815                  "'breakpoint_rule' in Piecewise object '%s'" % tmpname
816        if self._slope_component is not None and \
817               self._slope_rule is None:
818            raise TypeError, "Cannot specify 'slope_sets' without " \
819                  "'slope_rule' in Piecewise object '%s'" % tmpname
820
821
822        # Have we added implicit SOS constraints?
823        self._constraints_added = False
824
825    def construct(self, *args):
826        """
827        Construct the breakpoint and slope expressions
828        """
829
830        if pyomo.debug("verbose"):
831           print "Constructing Piecewise component, name="+self.name+", from args="+`args`
832
833        # NOTE: the _breakvals and _slopevals attributes can be either lists (if the
834        #       component is not indexed) or list of lists (if the component is indexed).
835
836        # Construct the breakpoints
837        self._breakvals = list()
838        if self._breakpoints is not None:
839            # Unpack breakpoints from an iterable
840            for exp in self._breakpoints:
841                if exp is not None:
842                    self._breakvals.append(exp)
843        else:
844            # Get breakpoints from a rule
845            if (self._breakpoint_component is None) or (self._breakpoint_component._ndim == 0):
846                # Rule isn't indexed
847                exp = self._breakpoint_rule((self.model,))
848                if not isinstance(exp, (list, tuple)):
849                   raise TypeError, "A breakpoint rule is expected to return a list or a tuple - returned type was "+str(type(exp))
850                if exp is not None:
851                    self._breakvals = exp
852            else:
853                # Rule is indexed
854                for index in self._breakpoint_component._index:
855                    # Make index a list, and then append self.model
856                    if isinstance(index, (list, tuple)):
857                        index = list(index)
858                    else:
859                        index = [index]
860                    index.append(self.model)
861
862                    # Get the next breakpoint
863                    exp = self._breakpoint_rule(*index)
864                    if exp is not None:
865                        self._breakvals.append(exp)
866
867        # Construct the slopes
868        self._slopevals = list()
869        if self._slopes is not None:
870            # Unpack slopes from an iterable
871            for exp in self._slopes:
872                if exp is not None:
873                    self._slopevals.append(exp)
874        else:
875            # Get slopes from a rule
876            if (self._slope_component is None) or (self._slope_component._ndim == 0):
877                # Rule isn't indexed
878                exp = self._slope_rule((self.model,))
879                if not isinstance(exp, (list, tuple)):
880                   raise TypeError, "A slope rule is expected to return a list or a tuple - returned type was "+str(type(exp))
881                if exp is not None:
882                    self._slopevals = exp
883            else:
884                # Rule is indexed
885                for index in self._slope_component._index:
886                    # Make index a list, and then append self.model
887                    if isinstance(index, (list, tuple)):
888                        index = list(index)
889                    else:
890                        index = [index]
891                    index.append(self.model)
892
893                    # Get the next breakpoint
894                    exp = self._slope_rule(*index)
895                    if exp is not None:
896                        self._slopevals.append(exp)
897
898        if pyomo.debug("verbose"):
899           if (self._slope_component is None) or (self._slope_component._ndim == 0):
900              print "Slopes="+str(self._slopevals)
901           else:
902              i = 0
903              for index in self._slope_component._index:
904                 print "Slope["+str(index)+"]="+str(self._slopevals[i])
905                 i += 1
906
907           if (self._breakpoint_component is None) or (self._breakpoint_component._ndim == 0):
908              print "Breakpoints="+str(self._slopevals)
909           else:
910              i = 0
911              for index in self._breakpoint_component._index:
912                 print "Breakpoint["+str(index)+"]="+str(self._breakvals[i])
913                 i += 1
914
915        # If N breakpoints are specified, then there can be N-1 or N+1 slopes.
916        if (self._breakpoint_component is None) or (self._breakpoint_component._ndim == 0):
917           if abs(len(self._slopevals) - len(self._breakvals)) != 1:
918               msg = "Error constructing Piecewise object '%(name)s'; "          \
919                     "%(breaks)i breakpoints and %(slopes)i slopes were "        \
920                     "specified.\n\tGiven %(breaks)i breakpoints, there should " \
921                     "be %(breaks)i-1 or %(breaks)i+1 slopes."
922               raise ValueError, msg % {
923                 "breaks" : len(self._breakvals),
924                 "slopes" : len(self._slopevals),
925                 "name"   : self.name
926               }
927        else:
928           for i in range(0,len(self._breakvals)):
929              if abs(len(self._slopevals[i]) - len(self._breakvals[i])) != 1:
930                  msg = "Error constructing Piecewise object '%(name)s'; "          \
931                        "%(breaks)i breakpoints and %(slopes)i slopes were "        \
932                        "specified.\n\tGiven %(breaks)i breakpoints, there should " \
933                        "be %(breaks)i-1 or %(breaks)i+1 slopes."
934                  raise ValueError, msg % {
935                    "breaks" : len(self._breakvals[i]),
936                    "slopes" : len(self._slopevals[i]),
937                    "name"   : self.name
938                  }
939
940        # If a user specified an unbounded problem, add "None" to both ends of self._breakvals.
941        if (self._breakpoint_component is None) or (self._breakpoint_component._ndim == 0):
942           if len(self._breakvals) < len(self._slopevals):
943              self._breakvals.insert(0, None)
944              self._breakvals.append(None)
945        else:
946           for i in range(0,len(self._breakvals)):
947              if len(self._breakvals[i]) < len(self._slopevals[i]):
948                 self._breakvals[i].insert(0, None)
949                 self._breakvals[i].append(None)             
950
951        self._generate_constraints()
952
953    def _concave_or_convex(self):
954        """
955        Determines if the expressions forms a concave or convex expression.
956        Returns -1 if concave, 1 if convex, 2 if affine, and 0 otherwise.
957        """
958
959        # Must be convex if there are less than three segments
960        if len(self._slopevals) <= 2:
961            if len(self._slopevals) == 1:
962                # Affine
963                return 2
964            elif self._slopevals[0] < self._slopevals[1]:
965                # Convex
966                return 1
967            elif self._slopevals[0] == self._slopevals[1]:
968                # Affine
969                return 2
970            else:
971                # Concave
972                return -1
973
974        # Find first segment that doesn't have the same slope,
975        # get initial direction
976        start = 1
977        while start < len(self._slopevals):
978            if value(self._slopevals[start]) == \
979                   value(self._slopevals[start-1]):
980                start += 1
981            else:
982                direction = value(self._slopevals[start]) > \
983                            value(self._slopevals[start-1])
984                break
985
986        # If we didn't find a change in slope, we're affine
987        if start == len(self._slopevals) - 1:
988            # Affine
989            return 2
990
991        # Make sure the rest of the segments follow the same direction
992        for i in xrange(start, len(self._slopevals)):
993            if not (direction and value(self._slopevals[i]) >= \
994                    value(self._slopevals[i-1])):
995                return 0
996
997        if direction:
998            # Convex
999            return 1
1000        else:
1001            # Concave
1002            return -1
1003
1004    def is_convex(self):
1005        res = self._concave_or_convex()
1006        return res == 1 or res == 2
1007
1008    def is_concave(self):
1009        res = self._concave_or_convex()
1010        return res == -1 or res == 2
1011
1012    def is_affine(self):
1013        return self._concave_or_convex() == 2
1014
1015    def is_concave_or_convex(self):
1016        return self._concave_or_convex() != 0
1017
1018    def _generate_constraints(self):
1019        """
1020        Generate the constraints on the variable. If this is a convex
1021        (or concave) expression, then we simply minimize (or maximize)
1022        over a number of linear constraints. If not, we introduce SOS2
1023        variables to constrain the expression.
1024        """
1025
1026        # Error if imported at module level
1027        from constraint import Constraint, SOSConstraint
1028        from rangeset import RangeSet
1029
1030        if self._constraints_added:
1031            # We've already made the implicit constraints
1032            return
1033
1034        # The root of all implicit constraint, variable, and set names
1035        baseName = str(self.name) + "_implicit_"
1036
1037        # Make a partial closure that returns a rule
1038        # Arguments passed become closures
1039        def makeRule(slope, hOffset, vOffset):
1040            expr = slope*(self.variable - hOffset) + vOffset
1041            if convex:
1042                def rule(model):
1043                    return self >= expr
1044            else:
1045                def rule(model):
1046                    return self <= expr
1047            return rule
1048
1049        def makeBoundRule(value, greater=True):
1050            if greater:
1051                def rule(model):
1052                    return self.variable >= value
1053            else:
1054                def rule(model):
1055                    return self.variable <= value
1056            return rule
1057
1058        # Bound model.value if domain is bounded
1059        if (self._breakvals[0] is not None):
1060            self.model.__setattr__(
1061                baseName + "value_lower_bound",
1062                Constraint(rule=makeBoundRule(self._breakvals[0], True))
1063                )
1064        if (self._breakvals[-1] is not None):
1065            self.model.__setattr__(
1066                baseName + "value_upper_bound",
1067                Constraint(rule=makeBoundRule(self._breakvals[-1], False))
1068                )
1069
1070        convex = self.is_convex()
1071        concave = self.is_concave()
1072        if (convex or concave) and (self._forcesos2 is False):
1073            # Maximize (convex) or minimize (concave) over all
1074            # constraints We first find the interval containing zero
1075
1076            # Handle the case where the domain has no lower bound.
1077            if self._breakvals[0] is None:
1078                start = 1
1079
1080                # Line segment with a lower domain bound of -inf
1081                self.model.__setattr__(
1082                    baseName + "constraint[%i]" % 0,
1083                    Constraint(rule=makeRule(self._slopevals[0],
1084                                             self._breakvals[1],
1085                                             self.offset
1086                                             ))
1087                    )
1088            else:
1089                start = 0
1090
1091            # Create the rest of the constraints. Walk from the lower
1092            # to upper bounds of the domain, keeping track of the
1093            # vertical offset.
1094            offset = self.offset
1095            i = start
1096            stop = len(self._breakvals)-1
1097            if (self._breakvals[-1] is None):
1098                stop -= 1
1099
1100            while i <= stop:
1101                slope = self._slopevals[i]
1102                # Make constraint
1103                self.model.__setattr__(
1104                    baseName + "constraint[%i]" % i,
1105                    Constraint(rule=makeRule(slope,
1106                                             self._breakvals[i],
1107                                             offset))
1108                    )
1109
1110                # Update offset
1111                if i < stop:
1112                    offset += slope*(self._breakvals[i+1] - self._breakvals[i])
1113
1114                i += 1
1115        else:
1116            # Introduce SOS2 constraints
1117
1118            # IMPORTANT
1119            #
1120            # For now, we assume that the domain is bounded. It will
1121            # fail otherwise. Support for extended domains coming
1122            # soon.
1123
1124            # Hack to artificially bound the problem, if necessary
1125            if self._breakvals[0] is None:
1126                self._breakvals = self._breakvals[1:]
1127                self._slopevals = self._slopevals[1:]
1128            if self._breakvals[-1] is None:
1129                self._breakvals = self._breakvals[:-1]
1130                self._slopevals = self._slopevals[:-1]
1131
1132            # Walk the breakpoints of the function and find the associated
1133            # images
1134            images = [self.offset]
1135            i = 0
1136            stop = len(self._slopevals) - 1
1137            while i <= stop:
1138                images.append(
1139                    images[-1] + self._slopevals[i] * \
1140                    (self._breakvals[i+1] - self._breakvals[i])
1141                    )
1142                i += 1
1143
1144            # Make the new set indexing the variables
1145            setName = baseName + "SOS_index_set"
1146            self.model.__setattr__(
1147                setName,
1148                RangeSet(0, len(self._breakvals)-1)
1149                )
1150
1151            # Make the new variables
1152            varName = baseName + "SOS_variables"
1153            self.model.__setattr__(
1154                varName,
1155                Var(self.model.__getattribute__(setName))
1156                )
1157
1158            # Constrain variables
1159
1160            # It would be nice to just bound all the SOS variables in the
1161            # Var() class, but then if the variable doesn't appear anywhere
1162            # else in the LP file CPLEX fails, so we explicity constrain
1163            # them here
1164
1165            def makeBoundsRule(varName):
1166                def boundsRule(i, model):
1167                    return (0, model.__getattribute__(varName)[i], 1)
1168                return boundsRule
1169
1170            self.model.__setattr__(
1171                baseName + "bounds_constraint",
1172                Constraint(self.model.__getattribute__(setName),
1173                           rule=makeBoundsRule(varName))
1174                )
1175
1176            # The SOS2 variables sum to 1
1177            def makeSumRule(varName, setName):
1178                def sumRule(model):
1179                    sum(model.__getattribute__(varName)[i] for i in \
1180                        model.__getattribute__(setName)) == 1
1181                return sumRule
1182
1183            self.model.__setattr__(
1184                baseName + "sum_constraint",
1185                Constraint(rule=makeSumRule(varName, setName)))
1186
1187            # The variable self is a convex combination of the images
1188
1189            def makeImageConvexCombRule(varName, setName, images):
1190                def imageConvexCombRule(model):
1191                    return self == sum(
1192                        model.__getattribute__(varName)[i] * images[i] \
1193                        for i in model.__getattribute__(setName))
1194
1195                return imageConvexCombRule
1196
1197            self.model.__setattr__(
1198                baseName + "image_convex_combination_constraint",
1199                Constraint(rule=makeImageConvexCombRule(varName,
1200                                                        setName,
1201                                                        images))
1202                )
1203
1204            # The variable model.value is a convex combination of the
1205            # breakpoints
1206
1207            def makeDomainConvexCombRule(varname, setName, breakpoints):
1208                def domainConvexCombRule(model):
1209                    return self.variable == sum(
1210                        model.__getattribute__(varName)[i] * breakpoints[i] \
1211                        for i in model.__getattribute__(setName))
1212
1213                return domainConvexCombRule
1214
1215            self.model.__setattr__(
1216                baseName + "domain_convex_combination_constraint",
1217                Constraint(rule=makeDomainConvexCombRule(varName,
1218                                                         setName,
1219                                                         self._breakvals
1220                                                         )))
1221
1222            # Finally, nonzero auxialliary variables must be adjacent
1223            self.model.__setattr__(
1224                baseName + "SOS2_constraint",
1225                SOSConstraint(var=self.model.__getattribute__(varName),
1226                              set=self.model.__getattribute__(setName),
1227                              sos=2))
1228
1229        # Don't make the constraints again
1230        self._constraints_added = True
Note: See TracBrowser for help on using the repository browser.