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
11import pyomo.common.unittest as unittest
12
13from pyomo.opt import (
14    TerminationCondition, SolutionStatus, check_available_solvers,
15)
16import pyomo.environ as pyo
17import pyomo.kernel as pmo
18import sys
19
20diff_tol = 1e-3
21
22mosek_available = check_available_solvers('mosek_direct')
23
24@unittest.skipIf(not mosek_available ,
25                 "MOSEK's python bindings are not available")
26class MOSEKDirectTests(unittest.TestCase):
27
28    def setUp(self):
29        self.stderr = sys.stderr
30        sys.stderr = None
31
32    def tearDown(self):
33        sys.stderr = self.stderr
34
35    def test_interface_call(self):
36
37        interface_instance = type(pyo.SolverFactory('mosek_direct'))
38        alt_1 = pyo.SolverFactory('mosek')
39        alt_2 = pyo.SolverFactory('mosek', solver_io='python')
40        alt_3 = pyo.SolverFactory('mosek', solver_io='direct')
41        self.assertIsInstance(alt_1, interface_instance)
42        self.assertIsInstance(alt_2, interface_instance)
43        self.assertIsInstance(alt_3, interface_instance)
44
45    def test_infeasible_lp(self):
46
47        model = pyo.ConcreteModel()
48        model.X = pyo.Var(within=pyo.NonNegativeReals)
49        model.C1 = pyo.Constraint(expr=model.X == 1)
50        model.C2 = pyo.Constraint(expr=model.X == 2)
51        model.O = pyo.Objective(expr=model.X)
52
53        opt = pyo.SolverFactory("mosek_direct")
54        results = opt.solve(model)
55
56        self.assertIn(results.solver.termination_condition,
57                      (TerminationCondition.infeasible,
58                       TerminationCondition.infeasibleOrUnbounded))
59
60    def test_unbounded_lp(self):
61
62        model = pyo.ConcreteModel()
63        model.X = pyo.Var()
64        model.O = pyo.Objective(expr=model.X)
65
66        opt = pyo.SolverFactory("mosek_direct")
67        results = opt.solve(model)
68
69        self.assertIn(results.solver.termination_condition,
70                      (TerminationCondition.unbounded,
71                       TerminationCondition.infeasibleOrUnbounded))
72
73    def test_optimal_lp(self):
74
75        model = pyo.ConcreteModel()
76        model.X = pyo.Var(within=pyo.NonNegativeReals)
77        model.O = pyo.Objective(expr=model.X)
78
79        opt = pyo.SolverFactory("mosek_direct")
80        results = opt.solve(model, load_solutions=False)
81
82        self.assertEqual(results.solution.status,
83                         SolutionStatus.optimal)
84
85    def test_get_duals_lp(self):
86
87        model = pyo.ConcreteModel()
88        model.X = pyo.Var(within=pyo.NonNegativeReals)
89        model.Y = pyo.Var(within=pyo.NonNegativeReals)
90
91        model.C1 = pyo.Constraint(expr=2*model.X + model.Y >= 8)
92        model.C2 = pyo.Constraint(expr=model.X + 3*model.Y >= 6)
93
94        model.O = pyo.Objective(expr=model.X + model.Y)
95
96        opt = pyo.SolverFactory("mosek_direct")
97        results = opt.solve(model, suffixes=['dual'], load_solutions=False)
98
99        model.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)
100        model.solutions.load_from(results)
101
102        self.assertAlmostEqual(model.dual[model.C1], 0.4, 4)
103        self.assertAlmostEqual(model.dual[model.C2], 0.2, 4)
104
105    def test_infeasible_mip(self):
106
107        model = pyo.ConcreteModel()
108        model.X = pyo.Var(within=pyo.NonNegativeIntegers)
109        model.C1 = pyo.Constraint(expr=model.X == 1)
110        model.C2 = pyo.Constraint(expr=model.X == 2)
111        model.O = pyo.Objective(expr=model.X)
112
113        opt = pyo.SolverFactory("mosek_direct")
114        results = opt.solve(model)
115
116        self.assertIn(results.solver.termination_condition,
117                      (TerminationCondition.infeasibleOrUnbounded,
118                       TerminationCondition.infeasible))
119
120    def test_unbounded_mip(self):
121
122        model = pyo.AbstractModel()
123        model.X = pyo.Var(within=pyo.Integers)
124        model.O = pyo.Objective(expr=model.X)
125
126        instance = model.create_instance()
127        opt = pyo.SolverFactory("mosek_direct")
128        results = opt.solve(instance)
129
130        self.assertIn(results.solver.termination_condition,
131                      (TerminationCondition.unbounded,
132                       TerminationCondition.infeasibleOrUnbounded))
133
134    def test_optimal_mip(self):
135
136        model = pyo.ConcreteModel()
137        model.X = pyo.Var(within=pyo.NonNegativeIntegers)
138        model.O = pyo.Objective(expr=model.X)
139
140        opt = pyo.SolverFactory("mosek_direct")
141        results = opt.solve(model, load_solutions=False)
142
143        self.assertEqual(results.solution.status,
144                         SolutionStatus.optimal)
145
146    def test_conic(self):
147
148        model = pmo.block()
149        model.o = pmo.objective(0.0)
150        model.c = pmo.constraint(body=0.0,
151                                 rhs=1)
152
153        b = model.quadratic = pmo.block()
154        b.x = pmo.variable_tuple((pmo.variable(),
155                                  pmo.variable()))
156        b.r = pmo.variable(lb=0)
157        b.c = pmo.conic.quadratic(x=b.x,
158                                  r=b.r)
159        model.o.expr += b.r
160        model.c.body += b.r
161        del b
162
163        b = model.rotated_quadratic = pmo.block()
164        b.x = pmo.variable_tuple((pmo.variable(),
165                                  pmo.variable()))
166        b.r1 = pmo.variable(lb=0)
167        b.r2 = pmo.variable(lb=0)
168        b.c = pmo.conic.rotated_quadratic(x=b.x,
169                                          r1=b.r1,
170                                          r2=b.r2)
171        model.o.expr += b.r1 + b.r2
172        model.c.body += b.r1 + b.r2
173        del b
174
175        import mosek
176        if mosek.Env().getversion() >= (9, 0, 0):
177            b = model.primal_exponential = pmo.block()
178            b.x1 = pmo.variable(lb=0)
179            b.x2 = pmo.variable()
180            b.r = pmo.variable(lb=0)
181            b.c = pmo.conic.primal_exponential(x1=b.x1,
182                                               x2=b.x2,
183                                               r=b.r)
184            model.o.expr += b.r
185            model.c.body += b.r
186            del b
187
188            b = model.primal_power = pmo.block()
189            b.x = pmo.variable_tuple((pmo.variable(),
190                                      pmo.variable()))
191            b.r1 = pmo.variable(lb=0)
192            b.r2 = pmo.variable(lb=0)
193            b.c = pmo.conic.primal_power(x=b.x,
194                                         r1=b.r1,
195                                         r2=b.r2,
196                                         alpha=0.6)
197            model.o.expr += b.r1 + b.r2
198            model.c.body += b.r1 + b.r2
199            del b
200
201            b = model.dual_exponential = pmo.block()
202            b.x1 = pmo.variable()
203            b.x2 = pmo.variable(ub=0)
204            b.r = pmo.variable(lb=0)
205            b.c = pmo.conic.dual_exponential(x1=b.x1,
206                                             x2=b.x2,
207                                             r=b.r)
208            model.o.expr += b.r
209            model.c.body += b.r
210            del b
211
212            b = model.dual_power = pmo.block()
213            b.x = pmo.variable_tuple((pmo.variable(),
214                                      pmo.variable()))
215            b.r1 = pmo.variable(lb=0)
216            b.r2 = pmo.variable(lb=0)
217            b.c = pmo.conic.dual_power(x=b.x,
218                                       r1=b.r1,
219                                       r2=b.r2,
220                                       alpha=0.4)
221            model.o.expr += b.r1 + b.r2
222            model.c.body += b.r1 + b.r2
223
224        opt = pmo.SolverFactory("mosek_direct")
225        results = opt.solve(model)
226
227        self.assertEqual(results.solution.status,
228                         SolutionStatus.optimal)
229
230
231if __name__ == "__main__":
232    unittest.main()
233