1#  ___________________________________________________________________________
2#
3#  Pyomo: Python Optimization Modeling Objects
4#  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC
5#  Under the terms of Contract DE-NA0003525 with National Technology and
6#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
7#  rights in this software.
8#  This software is distributed under the 3-clause BSD License.
9#  ___________________________________________________________________________
10
11"""Tests for the GDPopt solver plugin."""
12import logging
13from math import fabs
14from os.path import join, normpath
15
16from io import StringIO
17
18import pyomo.common.unittest as unittest
19from pyomo.common.log import LoggingIntercept
20from pyomo.common.collections import Bunch
21from pyomo.common.fileutils import import_file
22from pyomo.contrib.gdpopt.GDPopt import GDPoptSolver
23from pyomo.contrib.gdpopt.data_class import GDPoptSolveData
24from pyomo.contrib.gdpopt.mip_solve import solve_linear_GDP
25from pyomo.contrib.gdpopt.util import is_feasible, time_code
26from pyomo.environ import ConcreteModel, Objective, SolverFactory, Var, value, Integers, Block, Constraint, maximize
27from pyomo.gdp import Disjunct, Disjunction
28from pyomo.contrib.mcpp.pyomo_mcpp import mcpp_available
29from pyomo.opt import TerminationCondition
30
31from pyomo.common.fileutils import PYOMO_ROOT_DIR
32exdir = normpath(join(PYOMO_ROOT_DIR, 'examples', 'gdp'))
33
34mip_solver = 'glpk'
35nlp_solver = 'ipopt'
36global_nlp_solver = 'baron'
37global_nlp_solver_args = dict()
38minlp_solver = 'baron'
39LOA_solvers = (mip_solver, nlp_solver)
40GLOA_solvers = (mip_solver, global_nlp_solver, minlp_solver)
41LOA_solvers_available = all(SolverFactory(s).available() for s in LOA_solvers)
42GLOA_solvers_available = all(SolverFactory(s).available() for s in GLOA_solvers)
43license_available = SolverFactory(global_nlp_solver).license_is_valid() if GLOA_solvers_available else False
44
45
46class TestGDPoptUnit(unittest.TestCase):
47    """Real unit tests for GDPopt"""
48
49    @unittest.skipUnless(SolverFactory(mip_solver).available(), "MIP solver not available")
50    def test_solve_linear_GDP_unbounded(self):
51        m = ConcreteModel()
52        m.GDPopt_utils = Block()
53        m.x = Var(bounds=(-1, 10))
54        m.y = Var(bounds=(2, 3))
55        m.z = Var()
56        m.d = Disjunction(expr=[
57            [m.x + m.y >= 5], [m.x - m.y <= 3]
58        ])
59        m.o = Objective(expr=m.z)
60        m.GDPopt_utils.variable_list = [m.x, m.y, m.z]
61        m.GDPopt_utils.disjunct_list = [m.d._autodisjuncts[0], m.d._autodisjuncts[1]]
62        output = StringIO()
63        with LoggingIntercept(output, 'pyomo.contrib.gdpopt', logging.WARNING):
64            solver_data = GDPoptSolveData()
65            solver_data.timing = Bunch()
66            with time_code(solver_data.timing, 'main', is_main_timer=True):
67                solve_linear_GDP(m, solver_data, GDPoptSolver.CONFIG(dict(mip_solver=mip_solver, strategy='LOA')))
68            self.assertIn("Linear GDP was unbounded. Resolving with arbitrary bound values",
69                          output.getvalue().strip())
70
71    @unittest.skipUnless(SolverFactory(mip_solver).available(), "MIP solver not available")
72    def test_solve_lp(self):
73        m = ConcreteModel()
74        m.x = Var(bounds=(-5, 5))
75        m.c = Constraint(expr=m.x >= 1)
76        m.o = Objective(expr=m.x)
77        output = StringIO()
78        with LoggingIntercept(output, 'pyomo.contrib.gdpopt', logging.INFO):
79            SolverFactory('gdpopt').solve(m, mip_solver=mip_solver, strategy='LOA')
80            self.assertIn("Your model is an LP (linear program).",
81                          output.getvalue().strip())
82            self.assertAlmostEqual(value(m.o.expr), 1)
83
84    @unittest.skipUnless(SolverFactory(nlp_solver).available(), 'NLP solver not available')
85    def test_solve_nlp(self):
86        m = ConcreteModel()
87        m.x = Var(bounds=(-5, 5))
88        m.c = Constraint(expr=m.x >= 1)
89        m.o = Objective(expr=m.x ** 2)
90        output = StringIO()
91        with LoggingIntercept(output, 'pyomo.contrib.gdpopt', logging.INFO):
92            SolverFactory('gdpopt').solve(m, nlp_solver=nlp_solver, strategy='LOA')
93            self.assertIn("Your model is an NLP (nonlinear program).",
94                          output.getvalue().strip())
95            self.assertAlmostEqual(value(m.o.expr), 1)
96
97    @unittest.skipUnless(SolverFactory(mip_solver).available(), "MIP solver not available")
98    def test_solve_constant_obj(self):
99        m = ConcreteModel()
100        m.x = Var(bounds=(-5, 5))
101        m.c = Constraint(expr=m.x >= 1)
102        m.o = Objective(expr=1)
103        output = StringIO()
104        with LoggingIntercept(output, 'pyomo.contrib.gdpopt', logging.INFO):
105            SolverFactory('gdpopt').solve(m, mip_solver=mip_solver, strategy='LOA')
106            self.assertIn("Your model is an LP (linear program).",
107                          output.getvalue().strip())
108            self.assertAlmostEqual(value(m.o.expr), 1)
109
110    @unittest.skipUnless(SolverFactory(nlp_solver).available(), 'NLP solver not available')
111    def test_no_objective(self):
112        m = ConcreteModel()
113        m.x = Var(bounds=(-5, 5))
114        m.c = Constraint(expr=m.x ** 2 >= 1)
115        output = StringIO()
116        with LoggingIntercept(output, 'pyomo.contrib.gdpopt', logging.WARNING):
117            SolverFactory('gdpopt').solve(m, nlp_solver=nlp_solver, strategy='LOA')
118            self.assertIn("Model has no active objectives. Adding dummy objective.",
119                          output.getvalue().strip())
120
121    def test_multiple_objectives(self):
122        m = ConcreteModel()
123        m.x = Var()
124        m.o = Objective(expr=m.x)
125        m.o2 = Objective(expr=m.x + 1)
126        with self.assertRaisesRegex(ValueError, "Model has multiple active objectives"):
127            SolverFactory('gdpopt').solve(m, strategy='LOA')
128
129    def test_is_feasible_function(self):
130        m = ConcreteModel()
131        m.x = Var(bounds=(0, 3), initialize=2)
132        m.c = Constraint(expr=m.x == 2)
133        self.assertTrue(is_feasible(m, GDPoptSolver.CONFIG(dict(strategy='LOA'))))
134
135        m.c2 = Constraint(expr=m.x <= 1)
136        self.assertFalse(is_feasible(m, GDPoptSolver.CONFIG(dict(strategy='LOA'))))
137
138        m = ConcreteModel()
139        m.x = Var(bounds=(0, 3), initialize=2)
140        m.c = Constraint(expr=m.x >= 5)
141        self.assertFalse(is_feasible(m, GDPoptSolver.CONFIG(dict(strategy='LOA'))))
142
143        m = ConcreteModel()
144        m.x = Var(bounds=(3, 3), initialize=2)
145        self.assertFalse(is_feasible(m, GDPoptSolver.CONFIG(dict(strategy='LOA'))))
146
147        m = ConcreteModel()
148        m.x = Var(bounds=(0, 1), initialize=2)
149        self.assertFalse(is_feasible(m, GDPoptSolver.CONFIG(dict(strategy='LOA'))))
150
151        m = ConcreteModel()
152        m.x = Var(bounds=(0, 1), initialize=2)
153        m.d = Disjunct()
154        with self.assertRaisesRegex(NotImplementedError, "Found active disjunct"):
155            is_feasible(m, GDPoptSolver.CONFIG(dict(strategy='LOA')))
156
157
158@unittest.skipIf(not LOA_solvers_available,
159                 "Required subsolvers %s are not available"
160                 % (LOA_solvers,))
161class TestGDPopt(unittest.TestCase):
162    """Tests for the GDPopt solver plugin."""
163
164    def test_infeasible_GDP(self):
165        """Test for infeasible GDP."""
166        m = ConcreteModel()
167        m.x = Var(bounds=(0, 2))
168        m.d = Disjunction(expr=[
169            [m.x ** 2 >= 3, m.x >= 3],
170            [m.x ** 2 <= -1, m.x <= -1]])
171        m.o = Objective(expr=m.x)
172        output = StringIO()
173        with LoggingIntercept(output, 'pyomo.contrib.gdpopt', logging.WARNING):
174            SolverFactory('gdpopt').solve(
175                m, strategy='LOA',
176                mip_solver=mip_solver,
177                nlp_solver=nlp_solver)
178            self.assertIn("Set covering problem was infeasible.",
179                          output.getvalue().strip())
180
181    def test_GDP_nonlinear_objective(self):
182        m = ConcreteModel()
183        m.x = Var(bounds=(-1, 10))
184        m.y = Var(bounds=(2, 3))
185        m.d = Disjunction(expr=[
186            [m.x + m.y >= 5], [m.x - m.y <= 3]
187        ])
188        m.o = Objective(expr=m.x ** 2)
189        SolverFactory('gdpopt').solve(
190            m, strategy='LOA',
191            mip_solver=mip_solver,
192            nlp_solver=nlp_solver
193        )
194        self.assertAlmostEqual(value(m.o), 0)
195
196        m = ConcreteModel()
197        m.x = Var(bounds=(-1, 10))
198        m.y = Var(bounds=(2, 3))
199        m.d = Disjunction(expr=[
200            [m.x + m.y >= 5], [m.x - m.y <= 3]
201        ])
202        m.o = Objective(expr=-m.x ** 2, sense=maximize)
203        SolverFactory('gdpopt').solve(
204            m, strategy='LOA',
205            mip_solver=mip_solver,
206            nlp_solver=nlp_solver
207        )
208        self.assertAlmostEqual(value(m.o), 0)
209
210    def test_LOA_8PP_default_init(self):
211        """Test logic-based outer approximation with 8PP."""
212        exfile = import_file(
213            join(exdir, 'eight_process', 'eight_proc_model.py'))
214        eight_process = exfile.build_eight_process_flowsheet()
215        SolverFactory('gdpopt').solve(
216            eight_process, strategy='LOA',
217            mip_solver=mip_solver,
218            nlp_solver=nlp_solver,
219            tee=False)
220        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
221
222    @unittest.skipUnless(SolverFactory('gams').available(exception_flag=False), 'GAMS solver not available')
223    def test_LOA_8PP_gams_solver(self):
224        # Make sure that the duals are still correct
225        exfile = import_file(
226            join(exdir, 'eight_process', 'eight_proc_model.py'))
227        eight_process = exfile.build_eight_process_flowsheet()
228        SolverFactory('gdpopt').solve(
229            eight_process, strategy='LOA',
230            mip_solver=mip_solver,
231            nlp_solver='gams',
232            max_slack=0,
233            tee=False)
234        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
235
236    def test_LOA_8PP_force_NLP(self):
237        exfile = import_file(
238            join(exdir, 'eight_process', 'eight_proc_model.py'))
239        eight_process = exfile.build_eight_process_flowsheet()
240        SolverFactory('gdpopt').solve(
241            eight_process, strategy='LOA',
242            mip_solver=mip_solver,
243            nlp_solver=nlp_solver,
244            force_subproblem_nlp=True,
245            tee=False)
246        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
247
248    def test_LOA_strip_pack_default_init(self):
249        """Test logic-based outer approximation with strip packing."""
250        exfile = import_file(
251            join(exdir, 'strip_packing', 'strip_packing_concrete.py'))
252        strip_pack = exfile.build_rect_strip_packing_model()
253        SolverFactory('gdpopt').solve(
254            strip_pack, strategy='LOA',
255            mip_solver=mip_solver,
256            nlp_solver=nlp_solver)
257        self.assertTrue(
258            fabs(value(strip_pack.total_length.expr) - 11) <= 1E-2)
259
260    @unittest.category('expensive')
261    def test_LOA_constrained_layout_default_init(self):
262        """Test LOA with constrained layout."""
263        exfile = import_file(
264            join(exdir, 'constrained_layout', 'cons_layout_model.py'))
265        cons_layout = exfile.build_constrained_layout_model()
266        SolverFactory('gdpopt').solve(
267            cons_layout, strategy='LOA',
268            mip_solver=mip_solver,
269            nlp_solver=nlp_solver,
270            iterlim=120,
271            max_slack=5,  # problem is convex, so can decrease slack
272        )
273        objective_value = value(cons_layout.min_dist_cost.expr)
274        self.assertTrue(
275            fabs(objective_value - 41573) <= 200,
276            "Objective value of %s instead of 41573" % objective_value)
277
278    def test_LOA_8PP_maxBinary(self):
279        """Test logic-based OA with max_binary initialization."""
280        exfile = import_file(
281            join(exdir, 'eight_process', 'eight_proc_model.py'))
282        eight_process = exfile.build_eight_process_flowsheet()
283        SolverFactory('gdpopt').solve(
284            eight_process, strategy='LOA', init_strategy='max_binary',
285            mip_solver=mip_solver,
286            nlp_solver=nlp_solver)
287
288        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
289
290    def test_LOA_strip_pack_maxBinary(self):
291        """Test LOA with strip packing using max_binary initialization."""
292        exfile = import_file(
293            join(exdir, 'strip_packing', 'strip_packing_concrete.py'))
294        strip_pack = exfile.build_rect_strip_packing_model()
295        SolverFactory('gdpopt').solve(
296            strip_pack, strategy='LOA', init_strategy='max_binary',
297            mip_solver=mip_solver,
298            nlp_solver=nlp_solver)
299        self.assertTrue(
300            fabs(value(strip_pack.total_length.expr) - 11) <= 1E-2)
301
302    def test_LOA_8PP_fixed_disjuncts(self):
303        """Test LOA with 8PP using fixed disjuncts initialization."""
304        exfile = import_file(
305            join(exdir, 'eight_process', 'eight_proc_model.py'))
306        eight_process = exfile.build_eight_process_flowsheet()
307        initialize = [
308            # Use units 1, 4, 7, 8
309            eight_process.use_unit_1or2.disjuncts[0],
310            eight_process.use_unit_3ornot.disjuncts[1],
311            eight_process.use_unit_4or5ornot.disjuncts[0],
312            eight_process.use_unit_6or7ornot.disjuncts[1],
313            eight_process.use_unit_8ornot.disjuncts[0]
314        ]
315        for disj in eight_process.component_data_objects(Disjunct):
316            if disj in initialize:
317                disj.indicator_var.set_value(1)
318            else:
319                disj.indicator_var.set_value(0)
320        SolverFactory('gdpopt').solve(
321            eight_process, strategy='LOA', init_strategy='fix_disjuncts',
322            mip_solver=mip_solver,
323            nlp_solver=nlp_solver)
324
325        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
326
327    def test_LOA_custom_disjuncts(self):
328        """Test logic-based OA with custom disjuncts initialization."""
329        exfile = import_file(
330            join(exdir, 'eight_process', 'eight_proc_model.py'))
331        eight_process = exfile.build_eight_process_flowsheet()
332        initialize = [
333            # Use units 1, 4, 7, 8
334            [eight_process.use_unit_1or2.disjuncts[0],
335             eight_process.use_unit_3ornot.disjuncts[1],
336             eight_process.use_unit_4or5ornot.disjuncts[0],
337             eight_process.use_unit_6or7ornot.disjuncts[1],
338             eight_process.use_unit_8ornot.disjuncts[0]],
339            # Use units 2, 4, 6, 8
340            [eight_process.use_unit_1or2.disjuncts[1],
341             eight_process.use_unit_3ornot.disjuncts[1],
342             eight_process.use_unit_4or5ornot.disjuncts[0],
343             eight_process.use_unit_6or7ornot.disjuncts[0],
344             eight_process.use_unit_8ornot.disjuncts[0]]
345        ]
346
347        def assert_correct_disjuncts_active(nlp_model, solve_data):
348            if solve_data.master_iteration >= 1:
349                return  # only checking initialization
350            iter_num = solve_data.nlp_iteration
351            disjs_should_be_active = initialize[iter_num - 1]
352            for orig_disj, soln_disj in zip(
353                solve_data.original_model.GDPopt_utils.disjunct_list,
354                nlp_model.GDPopt_utils.disjunct_list
355            ):
356                if orig_disj in disjs_should_be_active:
357                    self.assertTrue(soln_disj.indicator_var.value == 1)
358
359        SolverFactory('gdpopt').solve(
360            eight_process, strategy='LOA', init_strategy='custom_disjuncts',
361            custom_init_disjuncts=initialize,
362            mip_solver=mip_solver,
363            nlp_solver=nlp_solver,
364            call_after_subproblem_feasible=assert_correct_disjuncts_active)
365
366        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
367
368
369@unittest.skipIf(not LOA_solvers_available,
370                 "Required subsolvers %s are not available"
371                 % (LOA_solvers,))
372class TestGDPoptRIC(unittest.TestCase):
373    """Tests for the GDPopt solver plugin."""
374
375    def test_infeasible_GDP(self):
376        """Test for infeasible GDP."""
377        m = ConcreteModel()
378        m.x = Var(bounds=(0, 2))
379        m.d = Disjunction(expr=[
380            [m.x ** 2 >= 3, m.x >= 3],
381            [m.x ** 2 <= -1, m.x <= -1]])
382        m.o = Objective(expr=m.x)
383        output = StringIO()
384        with LoggingIntercept(output, 'pyomo.contrib.gdpopt', logging.WARNING):
385            SolverFactory('gdpopt').solve(
386                m, strategy='RIC',
387                mip_solver=mip_solver,
388                nlp_solver=nlp_solver)
389            self.assertIn("Set covering problem was infeasible.",
390                          output.getvalue().strip())
391
392    def test_GDP_nonlinear_objective(self):
393        m = ConcreteModel()
394        m.x = Var(bounds=(-1, 10))
395        m.y = Var(bounds=(2, 3))
396        m.d = Disjunction(expr=[
397            [m.x + m.y >= 5], [m.x - m.y <= 3]
398        ])
399        m.o = Objective(expr=m.x ** 2)
400        SolverFactory('gdpopt').solve(
401            m, strategy='RIC',
402            mip_solver=mip_solver,
403            nlp_solver=nlp_solver
404        )
405        self.assertAlmostEqual(value(m.o), 0)
406
407        m = ConcreteModel()
408        m.x = Var(bounds=(-1, 10))
409        m.y = Var(bounds=(2, 3))
410        m.d = Disjunction(expr=[
411            [m.x + m.y >= 5], [m.x - m.y <= 3]
412        ])
413        m.o = Objective(expr=-m.x ** 2, sense=maximize)
414        SolverFactory('gdpopt').solve(
415            m, strategy='RIC',
416            mip_solver=mip_solver,
417            nlp_solver=nlp_solver
418        )
419        self.assertAlmostEqual(value(m.o), 0)
420
421    def test_RIC_8PP_default_init(self):
422        """Test logic-based outer approximation with 8PP."""
423        exfile = import_file(
424            join(exdir, 'eight_process', 'eight_proc_model.py'))
425        eight_process = exfile.build_eight_process_flowsheet()
426        SolverFactory('gdpopt').solve(
427            eight_process, strategy='RIC',
428            mip_solver=mip_solver,
429            nlp_solver=nlp_solver,
430            tee=False)
431        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
432
433    @unittest.skipUnless(SolverFactory('gams').available(exception_flag=False), 'GAMS solver not available')
434    def test_RIC_8PP_gams_solver(self):
435        # Make sure that the duals are still correct
436        exfile = import_file(
437            join(exdir, 'eight_process', 'eight_proc_model.py'))
438        eight_process = exfile.build_eight_process_flowsheet()
439        SolverFactory('gdpopt').solve(
440            eight_process, strategy='RIC',
441            mip_solver=mip_solver,
442            nlp_solver='gams',
443            max_slack=0,
444            tee=False)
445        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
446
447    def test_RIC_8PP_force_NLP(self):
448        exfile = import_file(
449            join(exdir, 'eight_process', 'eight_proc_model.py'))
450        eight_process = exfile.build_eight_process_flowsheet()
451        SolverFactory('gdpopt').solve(
452            eight_process, strategy='RIC',
453            mip_solver=mip_solver,
454            nlp_solver=nlp_solver,
455            force_subproblem_nlp=True,
456            tee=False)
457        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
458
459    def test_RIC_strip_pack_default_init(self):
460        """Test logic-based outer approximation with strip packing."""
461        exfile = import_file(
462            join(exdir, 'strip_packing', 'strip_packing_concrete.py'))
463        strip_pack = exfile.build_rect_strip_packing_model()
464        SolverFactory('gdpopt').solve(
465            strip_pack, strategy='RIC',
466            mip_solver=mip_solver,
467            nlp_solver=nlp_solver)
468        self.assertTrue(
469            fabs(value(strip_pack.total_length.expr) - 11) <= 1E-2)
470
471    @unittest.category('expensive')
472    def test_RIC_constrained_layout_default_init(self):
473        """Test RIC with constrained layout."""
474        exfile = import_file(
475            join(exdir, 'constrained_layout', 'cons_layout_model.py'))
476        cons_layout = exfile.build_constrained_layout_model()
477        SolverFactory('gdpopt').solve(
478            cons_layout, strategy='RIC',
479            mip_solver=mip_solver,
480            nlp_solver=nlp_solver,
481            iterlim=120,
482            max_slack=5,  # problem is convex, so can decrease slack
483        )
484        objective_value = value(cons_layout.min_dist_cost.expr)
485        self.assertTrue(
486            fabs(objective_value - 41573) <= 200,
487            "Objective value of %s instead of 41573" % objective_value)
488
489    def test_RIC_8PP_maxBinary(self):
490        """Test logic-based OA with max_binary initialization."""
491        exfile = import_file(
492            join(exdir, 'eight_process', 'eight_proc_model.py'))
493        eight_process = exfile.build_eight_process_flowsheet()
494        SolverFactory('gdpopt').solve(
495            eight_process, strategy='RIC', init_strategy='max_binary',
496            mip_solver=mip_solver,
497            nlp_solver=nlp_solver)
498
499        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
500
501    def test_RIC_strip_pack_maxBinary(self):
502        """Test RIC with strip packing using max_binary initialization."""
503        exfile = import_file(
504            join(exdir, 'strip_packing', 'strip_packing_concrete.py'))
505        strip_pack = exfile.build_rect_strip_packing_model()
506        SolverFactory('gdpopt').solve(
507            strip_pack, strategy='RIC', init_strategy='max_binary',
508            mip_solver=mip_solver,
509            nlp_solver=nlp_solver)
510        self.assertTrue(
511            fabs(value(strip_pack.total_length.expr) - 11) <= 1E-2)
512
513    def test_RIC_8PP_fixed_disjuncts(self):
514        """Test RIC with 8PP using fixed disjuncts initialization."""
515        exfile = import_file(
516            join(exdir, 'eight_process', 'eight_proc_model.py'))
517        eight_process = exfile.build_eight_process_flowsheet()
518        initialize = [
519            # Use units 1, 4, 7, 8
520            eight_process.use_unit_1or2.disjuncts[0],
521            eight_process.use_unit_3ornot.disjuncts[1],
522            eight_process.use_unit_4or5ornot.disjuncts[0],
523            eight_process.use_unit_6or7ornot.disjuncts[1],
524            eight_process.use_unit_8ornot.disjuncts[0]
525        ]
526        for disj in eight_process.component_data_objects(Disjunct):
527            if disj in initialize:
528                disj.indicator_var.set_value(1)
529            else:
530                disj.indicator_var.set_value(0)
531        SolverFactory('gdpopt').solve(
532            eight_process, strategy='RIC', init_strategy='fix_disjuncts',
533            mip_solver=mip_solver,
534            nlp_solver=nlp_solver)
535
536        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
537
538    def test_RIC_custom_disjuncts(self):
539        """Test logic-based OA with custom disjuncts initialization."""
540        exfile = import_file(
541            join(exdir, 'eight_process', 'eight_proc_model.py'))
542        eight_process = exfile.build_eight_process_flowsheet()
543        initialize = [
544            # Use units 1, 4, 7, 8
545            [eight_process.use_unit_1or2.disjuncts[0],
546             eight_process.use_unit_3ornot.disjuncts[1],
547             eight_process.use_unit_4or5ornot.disjuncts[0],
548             eight_process.use_unit_6or7ornot.disjuncts[1],
549             eight_process.use_unit_8ornot.disjuncts[0]],
550            # Use units 2, 4, 6, 8
551            [eight_process.use_unit_1or2.disjuncts[1],
552             eight_process.use_unit_3ornot.disjuncts[1],
553             eight_process.use_unit_4or5ornot.disjuncts[0],
554             eight_process.use_unit_6or7ornot.disjuncts[0],
555             eight_process.use_unit_8ornot.disjuncts[0]]
556        ]
557
558        def assert_correct_disjuncts_active(nlp_model, solve_data):
559            if solve_data.master_iteration >= 1:
560                return  # only checking initialization
561            iter_num = solve_data.nlp_iteration
562            disjs_should_be_active = initialize[iter_num - 1]
563            for orig_disj, soln_disj in zip(
564                solve_data.original_model.GDPopt_utils.disjunct_list,
565                nlp_model.GDPopt_utils.disjunct_list
566            ):
567                if orig_disj in disjs_should_be_active:
568                    self.assertTrue(soln_disj.indicator_var.value == 1)
569
570        SolverFactory('gdpopt').solve(
571            eight_process, strategy='RIC', init_strategy='custom_disjuncts',
572            custom_init_disjuncts=initialize,
573            mip_solver=mip_solver,
574            nlp_solver=nlp_solver,
575            call_after_subproblem_feasible=assert_correct_disjuncts_active)
576
577        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
578
579
580@unittest.skipIf(not GLOA_solvers_available,
581                 "Required subsolvers %s are not available"
582                 % (GLOA_solvers,))
583@unittest.skipIf(not mcpp_available(), "MC++ is not available")
584class TestGLOA(unittest.TestCase):
585    """Tests for global logic-based outer approximation."""
586
587    def test_GDP_integer_vars(self):
588        m = ConcreteModel()
589        m.x = Var(bounds=(0, 10))
590        m.y = Var(domain=Integers, bounds=(0, 5))
591        m.d = Disjunction(expr=[[
592            m.x >= m.y, m.y >= 3.5
593        ],
594        [
595            m.x >= m.y, m.y >= 2.5
596        ]])
597        m.o = Objective(expr=m.x)
598        SolverFactory('gdpopt').solve(
599            m, strategy='GLOA',
600            mip_solver=mip_solver,
601            nlp_solver=nlp_solver,
602            minlp_solver=minlp_solver
603        )
604        self.assertAlmostEqual(value(m.o.expr), 3)
605
606    def test_GDP_integer_vars_infeasible(self):
607        m = ConcreteModel()
608        m.x = Var(bounds=(0, 1))
609        m.y = Var(domain=Integers, bounds=(0, 5))
610        m.d = Disjunction(expr=[[
611            m.x >= m.y, m.y >= 3.5
612        ],
613        [
614            m.x >= m.y, m.y >= 2.5
615        ]])
616        m.o = Objective(expr=m.x)
617        res = SolverFactory('gdpopt').solve(
618            m, strategy='GLOA',
619            mip_solver=mip_solver,
620            nlp_solver=nlp_solver,
621            minlp_solver=minlp_solver
622        )
623        self.assertEqual(res.solver.termination_condition, TerminationCondition.infeasible)
624
625    @unittest.skipUnless(license_available, "Global NLP solver license not available.")
626    def test_GLOA_8PP(self):
627        """Test the global logic-based outer approximation algorithm."""
628        exfile = import_file(
629            join(exdir, 'eight_process', 'eight_proc_model.py'))
630        eight_process = exfile.build_eight_process_flowsheet()
631        SolverFactory('gdpopt').solve(
632            eight_process, strategy='GLOA', tee=False,
633            mip_solver=mip_solver,
634            nlp_solver=global_nlp_solver,
635            nlp_solver_args=global_nlp_solver_args
636        )
637        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
638
639    @unittest.skipUnless(license_available, "Global NLP solver license not available.")
640    def test_GLOA_8PP_force_NLP(self):
641        """Test the global logic-based outer approximation algorithm."""
642        exfile = import_file(
643            join(exdir, 'eight_process', 'eight_proc_model.py'))
644        eight_process = exfile.build_eight_process_flowsheet()
645        SolverFactory('gdpopt').solve(
646            eight_process, strategy='GLOA', tee=False,
647            mip_solver=mip_solver,
648            nlp_solver=global_nlp_solver,
649            nlp_solver_args=global_nlp_solver_args,
650            force_subproblem_nlp=True
651        )
652        self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1E-2)
653
654    @unittest.skipUnless(license_available, "Global NLP solver license not available.")
655    def test_GLOA_strip_pack_default_init(self):
656        """Test logic-based outer approximation with strip packing."""
657        exfile = import_file(
658            join(exdir, 'strip_packing', 'strip_packing_concrete.py'))
659        strip_pack = exfile.build_rect_strip_packing_model()
660        SolverFactory('gdpopt').solve(
661            strip_pack, strategy='GLOA',
662            mip_solver=mip_solver,
663            nlp_solver=global_nlp_solver,
664            nlp_solver_args=global_nlp_solver_args)
665        self.assertTrue(
666            fabs(value(strip_pack.total_length.expr) - 11) <= 1E-2)
667
668    @unittest.skipUnless(license_available, "Global NLP solver license not available.")
669    @unittest.category('expensive')
670    def test_GLOA_constrained_layout_default_init(self):
671        """Test LOA with constrained layout."""
672        exfile = import_file(
673            join(exdir, 'constrained_layout', 'cons_layout_model.py'))
674        cons_layout = exfile.build_constrained_layout_model()
675        SolverFactory('gdpopt').solve(
676            cons_layout, strategy='GLOA',
677            mip_solver=mip_solver,
678            iterlim=36,
679            nlp_solver=global_nlp_solver,
680            nlp_solver_args=global_nlp_solver_args,
681            tee=False)
682        objective_value = value(cons_layout.min_dist_cost.expr)
683        self.assertTrue(
684            fabs(objective_value - 41573) <= 200,
685            "Objective value of %s instead of 41573" % objective_value)
686
687    def test_GLOA_ex_633_trespalacios(self):
688        """Test LOA with Francisco thesis example."""
689        exfile = import_file(join(exdir, 'small_lit', 'ex_633_trespalacios.py'))
690        model = exfile.build_simple_nonconvex_gdp()
691        SolverFactory('gdpopt').solve(
692            model, strategy='GLOA',
693            mip_solver=mip_solver,
694            nlp_solver=global_nlp_solver,
695            nlp_solver_args=global_nlp_solver_args,
696            tee=False)
697        objective_value = value(model.obj.expr)
698        self.assertAlmostEqual(objective_value, 4.46, 2)
699
700    def test_GLOA_nonconvex_HENS(self):
701        exfile = import_file(join(exdir, 'small_lit', 'nonconvex_HEN.py'))
702        model = exfile.build_gdp_model()
703        SolverFactory('gdpopt').solve(
704            model, strategy='GLOA',
705            mip_solver=mip_solver,
706            nlp_solver=global_nlp_solver,
707            nlp_solver_args=global_nlp_solver_args,
708            tee=False)
709        objective_value = value(model.objective.expr)
710        self.assertAlmostEqual(objective_value * 1E-5, 1.14385, 2)
711
712    def test_GLOA_disjunctive_bounds(self):
713        exfile = import_file(join(exdir, 'small_lit', 'nonconvex_HEN.py'))
714        model = exfile.build_gdp_model()
715        SolverFactory('gdpopt').solve(
716            model, strategy='GLOA',
717            mip_solver=mip_solver,
718            nlp_solver=global_nlp_solver,
719            nlp_solver_args=global_nlp_solver_args,
720            calc_disjunctive_bounds=True,
721            tee=False)
722        objective_value = value(model.objective.expr)
723        self.assertAlmostEqual(objective_value * 1E-5, 1.14385, 2)
724
725
726if __name__ == '__main__':
727    unittest.main()
728