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
11from __future__ import division
12import pyomo.environ as pyo
13from pyomo.contrib.pynumero.interfaces.external_grey_box import \
14    ExternalGreyBoxBlock
15from pyomo.contrib.pynumero.examples.external_grey_box.react_example.reactor_model_outputs import ReactorConcentrationsOutputModel
16
17def maximize_cb_outputs(show_solver_log=False):
18    # in this simple example, we will use an external grey box model representing
19    # a steady-state reactor, and solve for the space velocity that maximizes
20    # the concentration of component B coming out of the reactor
21    m = pyo.ConcreteModel()
22
23    # create a block to store the external reactor model
24    m.reactor = ExternalGreyBoxBlock(
25        external_model=ReactorConcentrationsOutputModel()
26    )
27
28    # The reaction rate constants and the feed concentration will
29    # be fixed for this example
30    m.k1con = pyo.Constraint(expr=m.reactor.inputs['k1'] == 5/6)
31    m.k2con = pyo.Constraint(expr=m.reactor.inputs['k2'] == 5/3)
32    m.k3con = pyo.Constraint(expr=m.reactor.inputs['k3'] == 1/6000)
33    m.cafcon = pyo.Constraint(expr=m.reactor.inputs['caf'] == 10000)
34
35    # add an objective function that maximizes the concentration
36    # of cb coming out of the reactor
37    m.obj = pyo.Objective(expr=m.reactor.outputs['cb'], sense=pyo.maximize)
38
39    solver = pyo.SolverFactory('cyipopt')
40    solver.config.options['hessian_approximation'] = 'limited-memory'
41    results = solver.solve(m, tee=show_solver_log)
42    pyo.assert_optimal_termination(results)
43    return m
44
45if __name__ == '__main__':
46    m = maximize_cb_outputs(show_solver_log=True)
47    m.pprint()
48
49
50