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