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