1"""Functions for solving the master problem."""
2
3from __future__ import division
4
5from copy import deepcopy
6
7from pyomo.common.errors import InfeasibleConstraintException
8from pyomo.contrib.fbbt.fbbt import fbbt
9from pyomo.contrib.gdpopt.data_class import MasterProblemResult
10from pyomo.contrib.gdpopt.util import SuppressInfeasibleWarning, _DoNothing, get_main_elapsed_time
11from pyomo.core import (Block, Expression, Objective, TransformationFactory,
12                        Var, minimize, value, Constraint)
13from pyomo.gdp import Disjunct
14from pyomo.network import Port
15from pyomo.opt import SolutionStatus, SolverFactory
16from pyomo.opt import TerminationCondition as tc, SolverResults
17from pyomo.solvers.plugins.solvers.persistent_solver import PersistentSolver
18
19
20def solve_linear_GDP(linear_GDP_model, solve_data, config):
21    """Solves the linear GDP model and attempts to resolve solution issues."""
22    m = linear_GDP_model
23    GDPopt = m.GDPopt_utils
24    # Transform disjunctions
25    _bigm = TransformationFactory('gdp.bigm')
26    _bigm.handlers[Port] = False
27    _bigm.apply_to(m)
28
29    preprocessing_transformations = [
30        # # Propagate variable bounds
31        # 'contrib.propagate_eq_var_bounds',
32        # # Detect fixed variables
33        # 'contrib.detect_fixed_vars',
34        # # Propagate fixed variables
35        # 'contrib.propagate_fixed_vars',
36        # # Remove zero terms in linear expressions
37        # 'contrib.remove_zero_terms',
38        # # Remove terms in equal to zero summations
39        # 'contrib.propagate_zero_sum',
40        # # Transform bound constraints
41        # 'contrib.constraints_to_var_bounds',
42        # # Detect fixed variables
43        # 'contrib.detect_fixed_vars',
44        # # Remove terms in equal to zero summations
45        # 'contrib.propagate_zero_sum',
46        # Remove trivial constraints
47        'contrib.deactivate_trivial_constraints',
48    ]
49    if config.mip_presolve:
50        try:
51            fbbt(m, integer_tol=config.integer_tolerance)
52            for xfrm in preprocessing_transformations:
53                TransformationFactory(xfrm).apply_to(m)
54        except InfeasibleConstraintException:
55            config.logger.debug("MIP preprocessing detected infeasibility.")
56            mip_result = MasterProblemResult()
57            mip_result.feasible = False
58            mip_result.var_values = list(v.value for v in GDPopt.variable_list)
59            mip_result.pyomo_results = SolverResults()
60            mip_result.pyomo_results.solver.termination_condition = tc.error
61            mip_result.disjunct_values = list(
62                disj.binary_indicator_var.value
63                for disj in GDPopt.disjunct_list)
64            return mip_result
65
66    # Deactivate extraneous IMPORT/EXPORT suffixes
67    getattr(m, 'ipopt_zL_out', _DoNothing()).deactivate()
68    getattr(m, 'ipopt_zU_out', _DoNothing()).deactivate()
69
70    # Create solver, check availability
71    if not SolverFactory(config.mip_solver).available():
72        raise RuntimeError(
73            "MIP solver %s is not available." % config.mip_solver)
74
75    # Callback immediately before solving MIP master problem
76    config.call_before_master_solve(m, solve_data)
77
78    try:
79        with SuppressInfeasibleWarning():
80            mip_args = dict(config.mip_solver_args)
81            elapsed = get_main_elapsed_time(solve_data.timing)
82            remaining = max(config.time_limit - elapsed, 1)
83            if config.mip_solver == 'gams':
84                mip_args['add_options'] = mip_args.get('add_options', [])
85                mip_args['add_options'].append('option reslim=%s;' % remaining)
86            elif config.mip_solver == 'multisolve':
87                mip_args['time_limit'] = min(mip_args.get(
88                    'time_limit', float('inf')), remaining)
89            results = SolverFactory(config.mip_solver).solve(
90                m, **mip_args)
91    except RuntimeError as e:
92        if 'GAMS encountered an error during solve.' in str(e):
93            config.logger.warning(
94                "GAMS encountered an error in solve. Treating as infeasible.")
95            mip_result = MasterProblemResult()
96            mip_result.feasible = False
97            mip_result.var_values = list(v.value for v in GDPopt.variable_list)
98            mip_result.pyomo_results = SolverResults()
99            mip_result.pyomo_results.solver.termination_condition = tc.error
100            mip_result.disjunct_values = list(
101                disj.binary_indicator_var.value
102                for disj in GDPopt.disjunct_list)
103            return mip_result
104        else:
105            raise
106    terminate_cond = results.solver.termination_condition
107    if terminate_cond is tc.infeasibleOrUnbounded:
108        # Linear solvers will sometimes tell me that it's infeasible or
109        # unbounded during presolve, but fails to distinguish. We need to
110        # resolve with a solver option flag on.
111        results, terminate_cond = distinguish_mip_infeasible_or_unbounded(
112            m, config)
113    if terminate_cond is tc.unbounded:
114        # Solution is unbounded. Add an arbitrary bound to the objective and resolve.
115        # This occurs when the objective is nonlinear. The nonlinear objective is moved
116        # to the constraints, and deactivated for the linear master problem.
117        obj_bound = 1E15
118        config.logger.warning(
119            'Linear GDP was unbounded. '
120            'Resolving with arbitrary bound values of (-{0:.10g}, {0:.10g}) on the objective. '
121            'Check your initialization routine.'.format(obj_bound))
122        main_objective = next(m.component_data_objects(Objective, active=True))
123        GDPopt.objective_bound = Constraint(
124            expr=(-obj_bound, main_objective.expr, obj_bound))
125        with SuppressInfeasibleWarning():
126            results = SolverFactory(config.mip_solver).solve(
127                m, **config.mip_solver_args)
128        terminate_cond = results.solver.termination_condition
129
130    # Build and return results object
131    mip_result = MasterProblemResult()
132    mip_result.feasible = True
133    mip_result.var_values = list(v.value for v in GDPopt.variable_list)
134    mip_result.pyomo_results = results
135    mip_result.disjunct_values = list(
136        disj.binary_indicator_var.value for disj in GDPopt.disjunct_list)
137
138    if terminate_cond in {tc.optimal, tc.locallyOptimal, tc.feasible}:
139        pass
140    elif terminate_cond is tc.infeasible:
141        config.logger.info(
142            'Linear GDP is now infeasible. '
143            'GDPopt has finished exploring feasible discrete configurations.')
144        mip_result.feasible = False
145    elif terminate_cond is tc.maxTimeLimit:
146        # TODO check that status is actually ok and everything is feasible
147        config.logger.info(
148            'Unable to optimize linear GDP problem within time limit. '
149            'Using current solver feasible solution.')
150    elif (terminate_cond is tc.other and
151          results.solution.status is SolutionStatus.feasible):
152        # load the solution and suppress the warning message by setting
153        # solver status to ok.
154        config.logger.info(
155            'Linear GDP solver reported feasible solution, '
156            'but not guaranteed to be optimal.')
157    else:
158        raise ValueError(
159            'GDPopt unable to handle linear GDP '
160            'termination condition '
161            'of %s. Solver message: %s' %
162            (terminate_cond, results.solver.message))
163
164    return mip_result
165
166
167def distinguish_mip_infeasible_or_unbounded(m, config):
168    """Distinguish between an infeasible or unbounded solution.
169
170    Linear solvers will sometimes tell me that a problem is infeasible or
171    unbounded during presolve, but not distinguish between the two cases. We
172    address this by solving again with a solver option flag on.
173
174    """
175    tmp_args = deepcopy(config.mip_solver_args)
176    if config.mip_solver == 'gurobi':
177        # This solver option is specific to Gurobi.
178        tmp_args['options'] = tmp_args.get('options', {})
179        tmp_args['options']['DualReductions'] = 0
180    mipopt = SolverFactory(config.mip_solver)
181    if isinstance(mipopt, PersistentSolver):
182        mipopt.set_instance(m)
183    with SuppressInfeasibleWarning():
184        results = mipopt.solve(m, **tmp_args)
185    termination_condition = results.solver.termination_condition
186    return results, termination_condition
187
188
189def solve_LOA_master(solve_data, config):
190    """Solve the augmented lagrangean outer approximation master problem."""
191    m = solve_data.linear_GDP.clone()
192    GDPopt = m.GDPopt_utils
193    solve_data.mip_iteration += 1
194    main_objective = next(m.component_data_objects(Objective, active=True))
195
196    if solve_data.active_strategy == 'LOA':
197        # Set up augmented Lagrangean penalty objective
198        main_objective.deactivate()
199        sign_adjust = 1 if main_objective.sense == minimize else -1
200        GDPopt.OA_penalty_expr = Expression(
201            expr=sign_adjust * config.OA_penalty_factor *
202            sum(v for v in m.component_data_objects(
203                ctype=Var, descend_into=(Block, Disjunct))
204                if v.parent_component().local_name == 'GDPopt_OA_slacks'))
205        GDPopt.oa_obj = Objective(
206            expr=main_objective.expr + GDPopt.OA_penalty_expr,
207            sense=main_objective.sense)
208
209        obj_expr = GDPopt.oa_obj.expr
210        base_obj_expr = main_objective.expr
211    elif solve_data.active_strategy in {'GLOA', 'RIC'}:
212        obj_expr = base_obj_expr = main_objective.expr
213    else:
214        raise ValueError('Unrecognized strategy: ' + solve_data.active_strategy)
215
216    mip_result = solve_linear_GDP(m, solve_data, config)
217    if mip_result.feasible:
218        if main_objective.sense == minimize:
219            solve_data.LB = max(value(obj_expr), solve_data.LB)
220        else:
221            solve_data.UB = min(value(obj_expr), solve_data.UB)
222        solve_data.iteration_log[
223            (solve_data.master_iteration,
224             solve_data.mip_iteration,
225             solve_data.nlp_iteration)
226        ] = (
227            value(obj_expr),
228            value(base_obj_expr),
229            mip_result.var_values
230        )
231        config.logger.info(
232            'ITER {:d}.{:d}.{:d}-MIP: OBJ: {:.10g}  LB: {:.10g}  UB: {:.10g}'.format(
233                solve_data.master_iteration,
234                solve_data.mip_iteration,
235                solve_data.nlp_iteration,
236                value(obj_expr),
237                solve_data.LB, solve_data.UB))
238    else:
239        # Master problem was infeasible.
240        if solve_data.master_iteration == 1:
241            config.logger.warning(
242                'GDPopt initialization may have generated poor '
243                'quality cuts.')
244        # set optimistic bound to infinity
245        if main_objective.sense == minimize:
246            solve_data.LB = float('inf')
247        else:
248            solve_data.UB = float('-inf')
249    # Call the MILP post-solve callback
250    config.call_after_master_solve(m, solve_data)
251
252    return mip_result
253