source: coopr.fdt/trunk/coopr/fdt/FDT.py @ 2221

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

Updates due to changes in as_expr() syntax.

File size: 19.0 KB
Line 
1
2from coopr.opt import *
3import copy
4from pyutilib.misc import Container, Options
5import bidict
6from pyutilib.component.core import implements
7from coopr.pyomo import *
8import math
9
10class FDT(object):
11
12    # Declare that this is an IOptSolver plugin
13    implements(IOptSolver)
14
15    def initialize(self, options, instance):
16        """
17        Initialize the solver with an Options object
18        """
19        self.options = options
20        if self.options.always_prune is None:
21            self.options.always_prune = True
22        if self.options.timelimit == 0:
23            self.options.timelimit=None
24        if self.options.tol is None:
25            self.options.tol = 1e-8
26        if self.options.lpsolver is None:
27            self.options.lpsolver = 'cplex'
28        if self.options.smanager_type is None:
29            self.options.smanager_type = 'serial'
30        if self.options.max_pool_size is None:
31            self.options.max_pool_size = instance.statistics.number_of_binary_variables
32        #
33        # Create solver
34        #
35        self.opt = SolverFactory( options.lpsolver )
36        if self.opt is None:
37            raise ValueError, "Problem constructing solver `"+str(options.lpsolver)+"'"
38        self.opt.keepFiles=options.keep_files or options.log
39        #
40        # Create solver manager
41        #
42        self.solver_mngr = SolverManagerFactory( options.smanager_type )
43        if self.solver_mngr is None:
44            raise ValueError, "Problem constructing solver manager `"+str(options.smanager_type)+"'"
45        #
46        # Print FDT solver configuration
47        #
48        if options.debug:
49            print "Created solver FDT"
50            print options.__str__(indent='    ')
51
52    def solve(self, *args, **kwds):
53        """
54        Solve the specified problem and return a SolverResults object
55        """
56        # TODO: handling n
57        # TODO: complement objective costs for binary variables
58        # TODO: ensure that objective is minimization
59        instance = args[0]
60        if len(args) == 2:
61            options = args[1]
62        else:
63            options = Options()
64        for key in kwds:
65            options[key] = kwds[key]
66        self.initialize(options, instance)
67        print "Starting FDT heuristic"
68        #
69        # Verify that the problem can be solved
70        #
71        results = self.presolve(instance)
72        if not results is None:
73            return results
74        #
75        # Compute the LP relaxation of the model instance
76        # Terminate if the problem is infeasible or unbounded
77        #
78        [root_lp_value, root_lp_results] = self.process_root_lp(instance)
79        if root_lp_results.solver.termination_condition == TerminationCondition.infeasible:
80            root_lp_results.solver.termination_message = "Infeasible LP relaxation of FDT problem"
81            instance.store_info(root_lp_results)
82            return root_lp_results
83        if root_lp_results.solver.termination_condition == TerminationCondition.unbounded:
84            root_lp_results.solver.termination_message = "Unbounded LP relaxation of FDT problem"
85            instance.store_info(root_lp_results)
86            return root_lp_results
87        #
88        # TODO: terminate if all binary variables are integer in the root LP?
89        #
90        if self.options.debug:
91            instance.pprint()
92        #
93        # Generate a reduced instance, where variables with value less than
94        # a proscribed tolerance are fixed to zero.
95        #
96        reduced_instance, varmap = self.eliminate_zero_variables(instance, self.options.tol)
97        self.setup_branching_lp(reduced_instance)
98        #
99        # Setup the initial pool
100        #
101        rootlp = reduced_instance.solution()
102        omega = {}
103        for i in self.branching_lp.vars:
104            omega[i] = rootlp[i-1]
105        n = len(omega)
106        # TODO: distinguish between continuous vars in obj and binary vars
107        keys = list(self.branching_lp.binary_vars)
108        keys.sort()
109        pool = [Container(values=omega, fixed=None, integer=self.is_integer(omega, keys, self.options.tol), branch_vars=keys)]
110        if self.options.debug:
111            # reduced_instance.display()
112            print "OMEGA",omega
113            print "VARMAP",varmap
114        #
115        # Iterate until all locations are fixed
116        #
117        for i in keys:
118            tmppool=[]
119            did_work=False
120            for node in pool:
121                node.fixed = i
122                # TODO: also append to pool if the ith value is 0 or 1
123                if node.integer:
124                    tmppool.append(node)
125                    continue
126                branchinglp = self.generate_branching_lp(node, varmap, keys, reduced_instance)
127                results = self.solver_mngr.solve(branchinglp, opt=self.opt, keepFiles=True, tee=self.options.tee, timelimit=self.options.timelimit, options=" ".join(self.options.lpsolver_options))
128                if options.debug:
129                    results.write()
130                #
131                # TODO: process results better here...
132                # Q: is it possible to get an infeasible branching LP?  Unbounded?
133                #
134                if self.bad_results(results):
135                    print "Terminating due to error in results"
136                    return None
137                did_work |= self.process_results(node, results, tmppool)
138            print "Y POOL",did_work,tmppool
139            if did_work:
140                if (self.options.always_prune or len(pool) > n):
141                    pool = self.prune(tmppool, n, omega)
142                else:
143                    pool = tmppool
144            print "POOL after branching on variable",i
145            print pool
146        print "FINAL POOL"
147        print pool
148        # TODO: process final pool to eliminate things that aren't feasible
149        #   for the constraint matrix
150        # iterate for all 1's in the solution: try to drive the 1's to 0
151        #
152        for item in pool:
153            #BUG: reduced_instance.load(item.values)
154            for i in xsequence(len(item.values)):
155                reduced_instance.variable(i).value = item.values[i]
156            cost = reduced_instance.cost()
157            print "HERE", item.values, cost, "ratio=%s" % (root_lp_value / cost)
158        #
159        # Setup the results object
160        #
161        results = SolverResults()
162        instance.store_info(root_lp_results)
163        results.solver.status = SolverStatus.ok
164        results.solver.termination_condition = TerminationCondition.other
165        root_lp_results.solver.termination_message = "Processed all FDT iterations."
166        for item in pool:
167            soln = results.solution.add()
168            for i in self.branching_lp.vars:
169                reduced_instance.variable(i).value = item.values[i]
170                soln.variable[i] = item.values[i]
171            for obj in reduced_instance.objectives():
172                soln.objective[obj] = reduced_instance.objective(obj)()
173        #
174        return results
175
176    def presolve(self, instance):
177        """ Verify that the problem can be solved """
178        if instance.statistics.number_of_integer_variables == 0:
179            return None
180        results = SolverResults()
181        instance.store_info(results)
182        results.solver.status = SolverStatus.error
183        results.solver.message = "Cannot apply the FDT solver to a problem with general integer variables"
184        results.solver.termination_condition = TerminationCondition.error
185        return results
186
187    def process_root_lp(self, instance):
188        """
189        Solve for the LP relaxation for the non-continuous variables.
190        Fix zero valued variables
191        """
192        #
193        # Generate and solve model with relaxed integrality constraints
194        #
195        tmp = apply_transformation('relax_integrality', instance)
196        results = self.solver_mngr.solve(tmp, opt=self.opt, tee=self.options.tee, 
197                                                timelimit=self.options.timelimit, 
198                                                options=" ".join(self.options.lpsolver_options))
199        if results == None:
200            raise ValueError, "SolverManager returned None"
201        if results.solver.status != SolverStatus.ok:
202            return [None, results]
203        if results.solver.termination_condition == TerminationCondition.infeasible or \
204           results.solver.termination_condition == TerminationCondition.unbounded:
205            return [None, results]
206        if self.options.debug:
207            results.write()
208        print "Cost for initial root lp:", results.solution().objective.f.value
209        instance.load(results)
210        return [eval(results.solution().objective.f.value), results]
211
212    def eliminate_zero_variables(self, instance, tolerance):
213        #
214        # Find the nonzero variables, and limit the analysis to these
215        #
216        for var in instance.variables().values():
217            if var.is_binary() and math.fabs(var.value) < tolerance:
218                var.value = 0.0
219                var.fixed = True
220        reduced_instance = apply_transformation('eliminate_fixed_vars', instance)
221        varmap=bidict.bidict()
222        for i in xsequence(len(reduced_instance.variables())):
223            varmap[ reduced_instance.variable(i).id ] = reduced_instance.variable(i)._old_id
224        return reduced_instance, varmap
225
226    def bad_results(self, results):
227        return False
228
229    def generate_branching_lp(self, node, varmap, binary_vars, root):
230        #
231        # Add y_sum_bounds constraints until we reach the fix_var binary variable
232        #
233        self.branching_lp.y_sum_bounded.clear()
234        for i in binary_vars:
235            if i is node.fixed:
236                break
237            self.branching_lp.y_sum_bounded.add(i, self.branching_lp.y[i,0]+self.branching_lp.y[i,1] <= node.values[i])
238        self.branching_lp.y_sum_bounded.construct()
239        #
240        # Fix the values of y at the fixed variable
241        #
242        self.branching_lp.y_fixed.clear()
243        self.branching_lp.y_fixed.add(0, self.branching_lp.y[i,0] == 0)
244        # TODO: address dominant feasibility here
245        #self.branching_lp.y_fixed.add(1, self.branching_lp.y[i,1] <= self.branching_lp.l[1])
246        self.branching_lp.y_fixed.add(1, self.branching_lp.y[i,1] == self.branching_lp.l[1])
247        self.branching_lp.y_fixed.construct()
248        self.branching_lp.presolve()
249        #
250        # Print the final LP reformulation
251        #                   
252        if self.options.debug:
253            self.branching_lp.pprint()
254        #
255        return self.branching_lp
256
257    def is_integer(self, values, binary_vars, tol):
258        """
259        Returns True if all of the values in the list 'values' are within
260        the specified tolerance of being either 0 or 1.
261        """
262        for i in binary_vars:
263            if values[i] > tol and values[i] < 1.0-tol:
264                return False
265        return True
266
267    def setup_branching_lp(self, instance):
268        """
269        Setup the model for the branching LP.  Most of this model does not
270        change from one iteration to the next, so we set it up once.
271        """
272        #
273        model = Model()
274        #
275        model.vars = RangeSet(1,instance.nvariables())
276        #
277        def binary_vars_init(model):
278            return [i for i in model.vars if instance.variable(i).is_binary()]
279        model.binary_vars = Set(initialize=binary_vars_init)
280        #
281        model.cons = RangeSet(1,instance.nconstraints())
282        #
283        def range_cons_init(model):
284            """ Returns the set of indeces of range constraints """
285            #print "X",instance._con.keys(), instance.statistics.number_of_constraints, instance.nconstraints()
286            ans = []
287            for j in model.cons:
288                #print "Y",j
289                c = instance.constraint(j)
290                if c.lower is None or c.upper is None:
291                    continue
292                if math.fabs(c.lower - c.upper) >= self.options.tol:
293                    ans.append(j)
294            return ans
295        model.range_cons = Set(initialize=range_cons_init)
296
297        #
298        model.L = Set(initialize=set([0,1]))
299        #
300        model.l = Var(model.L, bounds=(0.0,1.0))
301        #
302        def Packdom_rule(model):
303            return model.l[0]+model.l[1]
304        model.Packdom = Objective(rule=Packdom_rule, sense=maximize)
305
306        #
307        model.y = Var(model.vars, model.L, bounds=(0.0,None))
308        #
309        def y_upper_rule(i,l,model):
310            return model.y[i,l] <= model.l[l]
311        model.y_upper = Constraint(model.vars, model.L, rule=y_upper_rule)
312
313        #
314        model.s = Var(model.range_cons, model.L)
315        #
316        model.C = Constraint(model.cons, model.L)
317        #
318        model.s_upper = Constraint(model.range_cons, model.L)
319        #
320        model.s_lower = Constraint(model.range_cons, model.L)
321
322        #
323        model.y_fixed = Constraint(model.L)
324        #
325        model.y_sum_bounded = Constraint(model.vars)
326
327        #
328        # Create an instance of this model
329        #
330        self.branching_lp = model.create()
331
332        #
333        # Setup dictionaries that are used to instantiate the expressions in C
334        #
335        vars0 = {}
336        vars1 = {}
337        for i in self.branching_lp.vars:
338            vars0[i] = self.branching_lp.y[i,0]
339            vars1[i] = self.branching_lp.y[i,1]
340        #
341        # Explicitly add constraints for C, s_lower and s_upper.
342        # These constraints are derived from the canonical expression of the
343        # constraints in 'instance'
344        #
345        for j in self.branching_lp.cons:
346            c = instance.constraint(j)
347            repn = generate_canonical_repn( c.body )
348            # TODO: error checking here?
349            if c.lower is None and c.upper is None:
350                # Ignore constraints with no bounds
351                continue
352            if c.lower is None or c.upper is None: 
353                if c.lower is None:
354                    # a'x <= b
355                    val = c.upper
356                    if 0 in repn:
357                        val -= repn[0][None]
358                        del repn[0]
359                    self.branching_lp.C.add((j,0), as_expr(repn, vars=vars0) <= val * self.branching_lp.l[0])
360                    self.branching_lp.C.add((j,1), as_expr(repn, vars=vars1) <= val * self.branching_lp.l[1])
361                else:
362                    # a'x >= b
363                    val = c.lower
364                    if 0 in repn:
365                        val -= repn[0][None]
366                        del repn[0]
367                    self.branching_lp.C.add((j,0), as_expr(repn, vars=vars0) >= val * self.branching_lp.l[0])
368                    self.branching_lp.C.add((j,1), as_expr(repn, vars=vars1) >= val * self.branching_lp.l[1])
369            else:
370                if math.fabs(c.lower - c.upper) < self.options.tol:
371                    # a'x = b
372                    val = c.lower
373                    if 0 in repn:
374                        val -= repn[0][None]
375                        del repn[0]
376                    self.branching_lp.C.add((j,0), as_expr(repn, vars=vars0) == val * self.branching_lp.l[0])
377                    self.branching_lp.C.add((j,1), as_expr(repn, vars=vars1) == val * self.branching_lp.l[1])
378                else:
379                    # b_l <= a'x <= b_u
380                    lower = c.lower
381                    upper = c.upper
382                    if 0 in repn:
383                        lower -= repn[0][None]
384                        upper -= repn[0][None]
385                        del repn[0]
386                    self.branching_lp.C.add((j,0), as_expr(repn, vars=vars0) == self.branching_lp.s[j,0])
387                    self.branching_lp.C.add((j,1), as_expr(repn, vars=vars1) == self.branching_lp.s[j,1])
388                    self.branching_lp.s_upper.add((j,0), self.branching_lp.s[j,0] <= upper * self.branching_lp.l[0])
389                    self.branching_lp.s_upper.add((j,1), self.branching_lp.s[j,1] <= upper * self.branching_lp.l[1])
390                    self.branching_lp.s_lower.add((j,0), self.branching_lp.s[j,0] >= lower * self.branching_lp.l[0])
391                    self.branching_lp.s_lower.add((j,1), self.branching_lp.s[j,1] >= lower * self.branching_lp.l[1])
392        self.branching_lp.C.construct()
393        self.branching_lp.s_upper.construct()
394        self.branching_lp.s_lower.construct()
395        #
396        # Print the final LP reformulation
397        #                   
398        if self.options.debug:
399            self.branching_lp.pprint()
400
401    def process_results(self, node, results, tmppool):
402        """ Process a branching LP results object and possibly update the pool """
403        # TODO: round parent values based on branching decisions (->1)
404        soln = results.solution()
405        did_work=False
406        if self.options.debug:
407            print soln.variable.keys()
408        if soln.variable.l[0].value > self.options.tol:
409            values = {}
410            for i in self.branching_lp.vars:
411                values[i] = soln.variable.y[i,0].value / soln.variable.l[0].value
412            tmppool += [ Container(fixed=None, values=values, integer=self.is_integer(values, node.branch_vars, self.options.tol), branch_vars=node.branch_vars) ]
413            did_work=True
414        if soln.variable.l[1].value > self.options.tol:
415            values = {}
416            for i in self.branching_lp.vars:
417                values[i] = soln.variable.y[i,1].value / soln.variable.l[1].value
418            did_work=True
419            tmppool += [ Container(fixed=None, values=values, integer=self.is_integer(values, node.branch_vars, self.options.tol), branch_vars=node.branch_vars) ]
420        print "X POOL",tmppool
421        return did_work
422       
423    def prune(self, pool, n, omega):
424        print "Pool", pool
425        #
426        # Create a packing LP model
427        #
428        model = Model()
429        #
430        model.V = RangeSet(1,n)
431        #
432        def xs_rule(i, model):
433            return omega[i]
434        model.xs = Param(model.V, initialize=xs_rule)
435        #
436        model.Good = RangeSet(1,len(pool))
437        #
438        def Sols_rule(i, u, model):
439            return pool[i-1].values[u]
440        model.Sols = Param(model.Good, model.V, initialize=Sols_rule)
441        #
442        model.lam = Var(model.Good, bounds=(0.0,1.0))
443        #
444        def Prunedom_rule(model):
445            return summation(model.lam)
446        model.Prunedom = Objective(rule=Prunedom_rule, sense=maximize)
447        #
448        def ConvComb_rule(u, model):
449            return sum([ model.lam[i]*model.Sols[i,u] for i in model.Good]) <= model.xs[u]
450        model.ConvComb = Constraint(model.V, rule=ConvComb_rule)
451        #
452        # Construct and solve this model
453        #
454        instance = model.create()
455        results = self.solver_mngr.solve(instance, opt=self.opt, tee=self.options.tee, timelimit=self.options.timelimit, options=" ".join(self.options.lpsolver_options))
456        #
457        # Create new pool using the nonzero lambdas
458        #
459        newpool = []
460        soln = results.solution()
461        print "Root LP", omega
462        print "Pool", pool
463        for i in soln.variable.lam:
464            print "Packing Info:", i, soln.variable.lam[i].value, pool[i-1].values
465            if soln.variable.lam[i].value > 0.0:
466                newpool.append( pool[i-1] )
467        return newpool
468
469
470# Register the solver with the name 'fdt'
471SolverRegistration('fdt', FDT)
472
Note: See TracBrowser for help on using the repository browser.