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.environ as pyo 12from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock 13from pyomo.contrib.pynumero.examples.external_grey_box.react_example.reactor_model_residuals import ReactorModel, ReactorModelWithHessian, \ 14 ReactorModelNoOutputs, ReactorModelScaled, create_pyomo_reactor_model 15 16def maximize_cb_ratio_residuals_with_output(show_solver_log=False, additional_options={}): 17 # in this simple example, we will use an external grey box model representing 18 # a steady-state reactor, and solve for the space velocity that maximizes 19 # the ratio of B to the other components coming out of the reactor 20 # This example illustrates the use of "equality constraints" or residuals 21 # in the external grey box example as well as outputs 22 m = pyo.ConcreteModel() 23 24 # create a block to store the external reactor model 25 m.reactor = ExternalGreyBoxBlock(external_model=ReactorModel()) 26 27 # The feed concentration will be fixed for this example 28 m.cafcon = pyo.Constraint(expr=m.reactor.inputs['caf'] == 10000) 29 30 # add an objective function that maximizes the concentration 31 # of cb coming out of the reactor 32 m.obj = pyo.Objective(expr=m.reactor.outputs['cb_ratio'], sense=pyo.maximize) 33 34 solver = pyo.SolverFactory('cyipopt') 35 solver.config.options['hessian_approximation'] = 'limited-memory' 36 for k,v in additional_options.items(): 37 solver.config.options[k] = v 38 results = solver.solve(m, tee=show_solver_log) 39 pyo.assert_optimal_termination(results) 40 return m 41 42def maximize_cb_ratio_residuals_with_hessian_with_output(show_solver_log=False, additional_options={}): 43 # in this simple example, we will use an external grey box model representing 44 # a steady-state reactor, and solve for the space velocity that maximizes 45 # the ratio of B to the other components coming out of the reactor 46 # This example illustrates the use of "equality constraints" or residuals 47 # in the external grey box example as well as outputs 48 m = pyo.ConcreteModel() 49 50 # create a block to store the external reactor model 51 m.reactor = ExternalGreyBoxBlock(external_model=ReactorModelWithHessian()) 52 53 # The feed concentration will be fixed for this example 54 m.cafcon = pyo.Constraint(expr=m.reactor.inputs['caf'] == 10000) 55 56 # add an objective function that maximizes the concentration 57 # of cb coming out of the reactor 58 m.obj = pyo.Objective(expr=m.reactor.outputs['cb_ratio'], sense=pyo.maximize) 59 60 solver = pyo.SolverFactory('cyipopt') 61 for k,v in additional_options.items(): 62 solver.config.options[k] = v 63 results = solver.solve(m, tee=show_solver_log) 64 pyo.assert_optimal_termination(results) 65 return m 66 67def maximize_cb_ratio_residuals_with_hessian_with_output_pyomo(show_solver_log=False, additional_options={}): 68 # this example is the same as the one above, but solves with a pure 69 # pyomo model - this is mostly for comparison and testing 70 m = create_pyomo_reactor_model() 71 solver = pyo.SolverFactory('ipopt') 72 for k,v in additional_options.items(): 73 solver.options[k] = v 74 solver.options['linear_solver'] = 'mumps' 75 results = solver.solve(m, tee=show_solver_log) 76 pyo.assert_optimal_termination(results) 77 return m 78 79def maximize_cb_ratio_residuals_with_output_scaling(show_solver_log=False, additional_options={}): 80 # in this simple example, we will use an external grey box model representing 81 # a steady-state reactor, and solve for the space velocity that maximizes 82 # the ratio of B to the other components coming out of the reactor 83 # This example illustrates the use of "equality constraints" or residuals 84 # in the external grey box example as well as outputs 85 86 # This example also shows how to do scaling. 87 # There are two things to scale, the "Pyomo" variables/constraints, 88 # and the external model residuals / output equations. 89 # The scaling factors for the external residuals and output equations 90 # are set in the derived ExternalGreyBoxModel class (see ReactorModelScaled 91 # for this example). 92 # The scaling factors for the Pyomo part of the model are set using suffixes 93 # - this requires that we declare the scaling suffix on the main model as 94 # shown below. 95 # - Then the scaling factors can be set directly on the Pyomo variables 96 # and constraints as also shown below. 97 # - Note: In this example, the scaling factors for the input and output 98 # variables from the grey box model are set in the finalize_block_construction 99 # callback (again, see ReactorModelScaled) 100 m = pyo.ConcreteModel() 101 102 # declare the scaling suffix on the model 103 m.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) 104 105 # create a block to store the external reactor model 106 m.reactor = ExternalGreyBoxBlock(external_model=ReactorModelScaled()) 107 108 # The feed concentration will be fixed for this example 109 m.cafcon = pyo.Constraint(expr=m.reactor.inputs['caf'] == 10000) 110 # set a scaling factor for this constraint - if we had additional pyomo 111 # variables, we could set them the same way 112 m.scaling_factor[m.cafcon] = 42.0 113 114 # add an objective function that maximizes the concentration 115 # of cb coming out of the reactor 116 m.obj = pyo.Objective(expr=m.reactor.outputs['cb_ratio'], sense=pyo.maximize) 117 118 solver = pyo.SolverFactory('cyipopt') 119 solver.config.options['hessian_approximation'] = 'limited-memory' 120 for k,v in additional_options.items(): 121 solver.config.options[k] = v 122 results = solver.solve(m, tee=show_solver_log) 123 pyo.assert_optimal_termination(results) 124 return m 125 126def maximize_cb_ratio_residuals_with_obj(show_solver_log=False, additional_options={}): 127 # in this simple example, we will use an external grey box model representing 128 # a steady-state reactor, and solve for the space velocity that maximizes 129 # the ratio of B to the other components coming out of the reactor 130 # This example illustrates the use of "equality constraints" or residuals 131 # in the external grey box example as well as additional pyomo variables 132 # and constraints 133 m = pyo.ConcreteModel() 134 135 # create a block to store the external reactor model 136 m.reactor = ExternalGreyBoxBlock() 137 m.reactor.set_external_model(ReactorModelNoOutputs()) 138 139 # The feed concentration will be fixed for this example 140 m.cafcon = pyo.Constraint(expr=m.reactor.inputs['caf'] == 10000) 141 142 # add an objective function that maximizes the concentration 143 # of cb coming out of the reactor 144 u = m.reactor.inputs 145 m.obj = pyo.Objective(expr=u['cb']/(u['ca']+u['cc']+u['cd']), sense=pyo.maximize) 146 147 solver = pyo.SolverFactory('cyipopt') 148 solver.config.options['hessian_approximation'] = 'limited-memory' 149 for k,v in additional_options.items(): 150 solver.config.options[k] = v 151 results = solver.solve(m, tee=show_solver_log) 152 pyo.assert_optimal_termination(results) 153 return m 154 155def maximize_cb_ratio_residuals_with_pyomo_variables(show_solver_log=False, additional_options={}): 156 # in this simple example, we will use an external grey box model representing 157 # a steady-state reactor, and solve for the space velocity that maximizes 158 # the ratio of B to the other components coming out of the reactor 159 # This example illustrates the use of "equality constraints" or residuals 160 # in the external grey box example as well as additional pyomo variables 161 # and constraints 162 m = pyo.ConcreteModel() 163 164 # create a block to store the external reactor model 165 m.reactor = ExternalGreyBoxBlock() 166 m.reactor.set_external_model(ReactorModelNoOutputs()) 167 168 # add a variable and constraint for the cb ratio 169 m.cb_ratio = pyo.Var(initialize=1) 170 u = m.reactor.inputs 171 m.cb_ratio_con = pyo.Constraint(expr = \ 172 u['cb']/(u['ca']+u['cc']+u['cd']) - m.cb_ratio == 0) 173 174 # The feed concentration will be fixed for this example 175 m.cafcon = pyo.Constraint(expr=m.reactor.inputs['caf'] == 10000) 176 177 # add an objective function that maximizes the concentration 178 # of cb coming out of the reactor 179 m.obj = pyo.Objective(expr=m.cb_ratio, sense=pyo.maximize) 180 181 solver = pyo.SolverFactory('cyipopt') 182 solver.config.options['hessian_approximation'] = 'limited-memory' 183 for k,v in additional_options.items(): 184 solver.config.options[k] = v 185 results = solver.solve(m, tee=show_solver_log) 186 pyo.assert_optimal_termination(results) 187 return m 188 189if __name__ == '__main__': 190 m = maximize_cb_ratio_residuals_with_output(show_solver_log=True) 191 192 # the next two are the same model with pyomo/ipopt and external/cyipopt 193 m = maximize_cb_ratio_residuals_with_hessian_with_output_pyomo(show_solver_log=True) 194 m = maximize_cb_ratio_residuals_with_hessian_with_output(show_solver_log=True) 195 196 aoptions={'hessian_approximation':'limited-memory', 197 #'limited_memory_update_type': 'sr1', 198 'nlp_scaling_method': 'user-scaling', 199 'print_level':10} 200 m = maximize_cb_ratio_residuals_with_output_scaling(show_solver_log=True, additional_options=aoptions) 201 202 m = maximize_cb_ratio_residuals_with_obj(show_solver_log=True) 203 204 m = maximize_cb_ratio_residuals_with_pyomo_variables(show_solver_log=True) 205