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#!/usr/bin/env python
12
13import pyomo.common.unittest as unittest
14from pyomo.common.config import ConfigBlock
15
16from pyomo.core.expr.current import identify_variables
17from pyomo.environ import ConcreteModel, Var, Reals, Objective, Constraint, ExternalFunction, SolverFactory, value, sqrt, sin
18from pyomo.opt import check_available_solvers
19
20try:
21    import numpy
22    numpy_available = True
23    from pyomo.contrib.trustregion.PyomoInterface import PyomoInterface
24except:
25    numpy_available = False
26
27
28@unittest.skipIf(numpy_available==False, "Skipping pyomo.contrib.trustregion tests because numpy is not installed.")
29class TestPyomoInterfaceInitialization(unittest.TestCase):
30    def setUp(self):
31        m = ConcreteModel()
32        m.z = Var(range(3),domain=Reals)
33        for i in range(3):
34            m.z[i] = 2.0
35        m.x = Var(range(2))
36        for i in range(2):
37            m.x[i] = 2.0
38        m.obj = Objective(expr= (m.z[0]-1.0)**2 + (m.z[0]-m.z[1])**2 + (m.z[2]-1.0)**2 + (m.x[0]-1.0)**4 + (m.x[1]-1.0)**6)
39        m.c2 = Constraint(expr=m.z[2]**4 * m.z[1]**2 + m.z[1] == 8+sqrt(2.0))
40        self.m = m
41
42    def test_1(self):
43        '''
44        The simplest case that the black box has only two inputs and there is only one black block involved
45        '''
46        def blackbox(a,b):
47            return sin(a-b)
48
49        m = self.m
50        bb = ExternalFunction(blackbox)
51        m.eflist = [bb]
52        m.c1 = Constraint(expr=m.x[0] * m.z[0]**2 + bb(m.x[0],m.x[1]) == 2*sqrt(2.0))
53        pI = PyomoInterface(m, [bb], ConfigBlock())
54        self.assertEqual(pI.lx,2)
55        self.assertEqual(pI.ly,1)
56        self.assertEqual(pI.lz,3)
57        self.assertEqual(len(list(identify_variables(m.c1.body))),3)
58        self.assertEqual(len(list(identify_variables(m.c2.body))),2)
59
60    def test_2(self):
61        '''
62        The simplest case that the black box has only one inputs and there is only a formula
63        '''
64        def blackbox(a):
65            return sin(a)
66
67        m = self.m
68        bb = ExternalFunction(blackbox)
69        m.eflist = [bb]
70        m.c1 = Constraint(expr=m.x[0] * m.z[0]**2 + bb(m.x[0]-m.x[1]) == 2*sqrt(2.0))
71        pI = PyomoInterface(m, [bb], ConfigBlock())
72        self.assertEqual(pI.lx,1)
73        self.assertEqual(pI.ly,1)
74        self.assertEqual(pI.lz,5)
75        self.assertEqual(len(list(identify_variables(m.c1.body))),3)
76        self.assertEqual(len(list(identify_variables(m.c2.body))),2)
77        self.assertEqual(len(m.tR.conset),1)
78        self.assertEqual(len(list(identify_variables(m.tR.conset[1].body))),3)
79
80    def test_3(self):
81        '''
82        The simplest case that the black box has only two inputs and there is only one black block involved
83        '''
84        def blackbox(a,b):
85            return sin(a-b)
86
87        m = self.m
88        bb = ExternalFunction(blackbox)
89        m.eflist = [bb]
90        m.c1 = Constraint(expr=m.x[0] * m.z[0]**2 + bb(m.x[0], m.x[1]) == 2*sqrt(2.0))
91        m.c3 = Constraint(expr=m.x[0] * m.z[0]**2 + bb(m.x[0], m.z[1]) == 2*sqrt(2.0))
92        pI = PyomoInterface(m, [bb], ConfigBlock())
93        self.assertEqual(pI.lx,3)
94        self.assertEqual(pI.ly,2)
95        self.assertEqual(pI.lz,2)
96        self.assertEqual(len(list(identify_variables(m.c1.body))),3)
97        self.assertEqual(len(list(identify_variables(m.c2.body))),2)
98        self.assertEqual(len(list(identify_variables(m.c3.body))),3)
99
100
101    def test_4(self):
102        '''
103        The simplest case that the black box has only two inputs and there is only one black block involved
104        '''
105        def blackbox(a,b):
106            return sin(a-b)
107
108        m = self.m
109        bb = ExternalFunction(blackbox)
110        m.eflist = [bb]
111        m.c1 = Constraint(expr=m.x[0] * m.z[0]**2 + bb(m.x[0], m.x[1]) == 2*sqrt(2.0))
112        m.c3 = Constraint(expr=m.x[0] * m.z[0]**2 + bb(m.x[0], m.z[1]) == 2*sqrt(2.0))
113        pI = PyomoInterface(m, [bb], ConfigBlock())
114        self.assertEqual(pI.lx,3)
115        self.assertEqual(pI.ly,2)
116        self.assertEqual(pI.lz,2)
117        self.assertEqual(len(list(identify_variables(m.c1.body))),3)
118        self.assertEqual(len(list(identify_variables(m.c2.body))),2)
119        self.assertEqual(len(list(identify_variables(m.c3.body))),3)
120
121    @unittest.skipIf(not check_available_solvers('ipopt'),
122                     "The 'ipopt' solver is not available")
123    @unittest.skipIf(not check_available_solvers('gjh'),
124                     "The 'gjh' solver is not available")
125    def test_execute_TRF(self):
126        m = ConcreteModel()
127        m.z = Var(range(3), domain=Reals, initialize=2.)
128        m.x = Var(range(2), initialize=2.)
129        m.x[1] = 1.0
130
131        def blackbox(a,b):
132            return sin(a-b)
133        bb = ExternalFunction(blackbox)
134
135        m.obj = Objective(
136            expr=(m.z[0]-1.0)**2 + (m.z[0]-m.z[1])**2 + (m.z[2]-1.0)**2 \
137            + (m.x[0]-1.0)**4 + (m.x[1]-1.0)**6 # + m.bb(m.x[0],m.x[1])
138        )
139        m.c1 = Constraint(
140            expr=m.x[0] * m.z[0]**2 + bb(m.x[0],m.x[1]) == 2*sqrt(2.0))
141        m.c2 = Constraint(expr=m.z[2]**4 * m.z[1]**2 + m.z[1] == 8+sqrt(2.0))
142
143        SolverFactory('trustregion').solve(m, [bb])
144
145        self.assertAlmostEqual(value(m.obj), 0.277044789315, places=4)
146        self.assertAlmostEqual(value(m.x[0]), 1.32193855369, places=4)
147        self.assertAlmostEqual(value(m.x[1]), 0.628744699822, places=4)
148
149
150if __name__ == '__main__':
151    unittest.main()
152