diff --git a/performance/models/implicit/__init__.py b/performance/models/implicit/__init__.py new file mode 100644 index 0000000..00f3bb9 --- /dev/null +++ b/performance/models/implicit/__init__.py @@ -0,0 +1,31 @@ +__doc__ = ( + """A directory of tests for the implicit function interface provided + by ExternalPyomoModel. + + Pyomo models interfaced with implicit functions via ExternalPyomoModel + are solved via a direct interface to a nonlinear solver, e.g. CyIpopt, + and there is not a single nl file that can encode the optimization + problem. We therefore track the following data during a performance + test: + + create_instance + setup_implicit + solve + check_results + + The "setup_implicit" label corresponds to the creation of a model with + ExternalGreyBoxBlock components containing the ExternalPyomoModel + objects that encode the implicit functions. + The "solve" label corresponds to the construction of a + PyomoNLPWithGreyBoxBlocks object and solution via CyIpopt with repeated + evaluation of the function and derivative methods of ExternalPyomoModel. + These methods involve heavy use of Pyomo components and expressions + (via calculate_variable_from_constraint), so including the solve should + be valuable for Pyomo performance testing. + It may be less valuable, however, if performance is found to be too + highly dependent on Ipopt/ASL versions. + Creating the ExternalPyomoModel and NLP objects also involve nl file + writes, which may be possible to track individually at a later date. + + """ +) diff --git a/performance/models/implicit/distill.dat b/performance/models/implicit/distill.dat new file mode 100644 index 0000000..c39307a --- /dev/null +++ b/performance/models/implicit/distill.dat @@ -0,0 +1,12 @@ +set S_TRAYS := +1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32; + +set S_RECTIFICATION := +2 3 4 5 6 7 8 9 10 11 12 13 14 15 16; + +set S_STRIPPING := +18 19 20 21 22 23 24 25 26 27 28 29 30 31; + +param x0 := +1 0.935419416 2 0.900525537 3 0.862296451 4 0.821699403 5 0.779990796 6 0.738571686 7 0.698804909 8 0.661842534 9 0.628507776 10 0.5992527 11 0.57418568 12 0.553144227 13 0.535784544 14 0.52166551 15 0.510314951 16 0.501275092 17 0.494128917 18 0.48544992 19 0.474202481 20 0.459803499 21 0.441642973 22 0.419191098 23 0.392055492 24 0.360245926 25 0.32407993 26 0.284676816 27 0.243209213 28 0.201815683 29 0.16177269 30 0.12514971 31 0.092458326 32 0.064583177; + diff --git a/performance/models/implicit/distill_DAE.py b/performance/models/implicit/distill_DAE.py new file mode 100644 index 0000000..08cf576 --- /dev/null +++ b/performance/models/implicit/distill_DAE.py @@ -0,0 +1,249 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import os +import pyomo.environ as pyo +from pyomo.environ import * +from pyomo.dae import * +from pyomo.util.subsystems import TemporarySubsystemManager +from pyomo.util.calc_var_value import calculate_variable_from_constraint +from pyomo.contrib.pynumero.interfaces.external_pyomo_model import ( + ExternalPyomoModel, +) +from pyomo.contrib.pynumero.interfaces.external_grey_box import ( + ExternalGreyBoxBlock, +) + +def make_model(horizon=52, vol=1.6, x_Feed=0.5): + + model = AbstractModel() + + model.Feed = Param(initialize = 24.0/60.0) + model.x_Feed = Param(initialize = x_Feed) + model.D = Param(initialize = model.x_Feed*model.Feed) + model.vol = Param(initialize = vol) + model.atray = Param(initialize = 0.25) + model.acond = Param(initialize = 0.5) + model.areb = Param(initialize = 1.0) + + model.S_TRAYS = Set(dimen=1) + model.S_RECTIFICATION = Set(within = model.S_TRAYS) + model.S_STRIPPING = Set(within = model.S_TRAYS) + model.x0 = Param(model.S_TRAYS) + + model.t = ContinuousSet(initialize=range(1, horizon)) + # Alternatively you could simply specify bounds on the + # ContinuousSet and let the finite element points be generated + # automatically. + # model.t = ContinuousSet(bounds=(1,51)) + + model.y = Var(model.S_TRAYS, model.t) + def x_init_rule(m,n,ti): + return value(m.x0[n]) + model.x = Var(model.S_TRAYS, model.t, initialize=x_init_rule) + model.dx = DerivativeVar(model.x) + + model.rr = Var(model.t,initialize=3.0) + model.L = Var(model.t,initialize=0.6) + model.V = Var(model.t,initialize=0.8) + model.FL = Var(model.t,initialize=1) + model.u1 = Var(model.t,initialize=3.0, bounds=(1,5)) + + model.alpha = Param(initialize = 1000) + model.rho = Param(initialize = 1) + model.u1_ref = Param(initialize = 2.0) + model.y1_ref = Param(initialize = 0.895814) + + ### + # Model constraints + ### + def reflux_ratio_rule(m,t): + return m.rr[t] == m.u1[t] + model.reflux_ratio = Constraint(model.t, rule=reflux_ratio_rule) + + def flowrate_rectification_rule(m,t): + return m.L[t] == m.rr[t]*m.D + model.flowrate_rectification = Constraint(model.t, rule=flowrate_rectification_rule) + + def vapor_column_rule(m,t): + return m.V[t] == m.L[t]+m.D + model.vapor_column = Constraint(model.t, rule=vapor_column_rule) + + def flowrate_stripping_rule(m,t): + return m.FL[t] == m.Feed + m.L[t] + model.flowrate_stripping = Constraint(model.t, rule=flowrate_stripping_rule) + + def mole_frac_balance_rule(m,n,t): + return m.y[n,t] == m.x[n,t]*m.vol/(1+((m.vol-1)*m.x[n,t])) + model.mole_frac_balance = Constraint(model.S_TRAYS, model.t, rule=mole_frac_balance_rule) + + def _diffeq(m,n,t): + + if t == 1: + return Constraint.Skip + if n == 1: + return m.dx[n,t] == 1/m.acond*m.V[t]*(m.y[n+1,t]-m.x[n,t]) + elif n in m.S_RECTIFICATION: + return m.dx[n,t] == 1/m.atray*(m.L[t]*(m.x[n-1,t]-m.x[n,t])-m.V[t]*(m.y[n,t]-m.y[n+1,t])) + elif n == 17: + return m.dx[n,t] == 1/m.atray*(m.Feed*m.x_Feed+m.L[t]*m.x[n-1,t]-m.FL[t]*m.x[n,t]-m.V[t]*(m.y[n,t]-m.y[n+1,t])) + elif n in m.S_STRIPPING: + return m.dx[n,t] == 1/m.atray*(m.FL[t]*(m.x[n-1,t]-m.x[n,t])-m.V[t]*(m.y[n,t]-m.y[n+1,t])) + else : + return m.dx[n,t] == 1/m.areb*(m.FL[t]*m.x[n-1,t]-(m.Feed-m.D)*m.x[n,t]-m.V[t]*m.y[n,t]) + model.diffeq = Constraint(model.S_TRAYS, model.t, rule=_diffeq) + + def _init_rule(m,n): + return m.x[n,1] == m.x0[n] + model.init_rule = Constraint(model.S_TRAYS, rule=_init_rule) + + return model + + +def discretize_model(instance): + # Discretize using Finite Difference Approach + discretizer = pyo.TransformationFactory('dae.finite_difference') + discretizer.apply_to(instance,nfe=50,scheme='BACKWARD') + + # Discretize using Orthogonal Collocation + # discretizer = TransformationFactory('dae.collocation') + # discretizer.apply_to(instance,nfe=50,ncp=3) + + # The objective function in the manually discretized pyomo model + # iterated over all finite elements and all collocation points. Since + # the objective function is not explicitly indexed by a ContinuousSet + # we add the objective function to the model after it has been + # discretized to ensure that we include all the discretization points + # when we take the sum. + + +def add_objective(instance): + def obj_rule(m): + return m.alpha*sum((m.y[1,i] - m.y1_ref)**2 for i in m.t if i != 1) + m.rho*sum((m.u1[i] - m.u1_ref)**2 for i in m.t if i!=1) + instance.OBJ = pyo.Objective(rule=obj_rule) + + # Calculate setpoint for the reduced space. + # Reduced space objective must not use algebraic variables. + t0 = instance.t.first() + to_reset = [instance.y[1, t0], instance.x[1, t0]] + with TemporarySubsystemManager(to_reset=to_reset): + instance.y[1, t0].set_value(pyo.value(instance.y1_ref)) + calculate_variable_from_constraint( + instance.x[1, t0], + instance.mole_frac_balance[1, t0], + ) + instance.x1_ref = pyo.Param(initialize=instance.x[1, t0].value) + + def rs_obj_rule(m): + return ( + m.alpha*sum((m.x[1,i] - m.x1_ref)**2 for i in m.t if i != 1) + + m.rho*sum((m.u1[i] - m.u1_ref)**2 for i in m.t if i != 1) + ) + instance.REDUCED_SPACE_OBJ = pyo.Objective(rule=rs_obj_rule) + instance.OBJ.deactivate() + + +def create_instance(): + model = make_model() + file_dir = os.path.dirname(__file__) + fname = os.path.join(file_dir, "distill.dat") + instance = model.create_instance(fname) + discretize_model(instance) + add_objective(instance) + return instance + + +def setup_implicit(instance): + diff_vars = [pyo.Reference(instance.x[i, :]) for i in instance.S_TRAYS] + deriv_vars = [pyo.Reference(instance.dx[i, :]) for i in instance.S_TRAYS] + disc_eqns = [pyo.Reference(instance.dx_disc_eq[i, :]) for i in instance.S_TRAYS] + diff_eqns = [pyo.Reference(instance.diffeq[i, :]) for i in instance.S_TRAYS] + + n_diff = len(diff_vars) + assert n_diff == len(deriv_vars) + assert n_diff == len(disc_eqns) + assert n_diff == len(diff_eqns) + + alg_vars = [] + alg_eqns = [] + alg_vars.extend(pyo.Reference(instance.y[i, :]) for i in instance.S_TRAYS) + alg_eqns.extend(pyo.Reference(instance.mole_frac_balance[i, :]) + for i in instance.S_TRAYS) + # Since we are not adding them to the reduced space model, alg vars do not + # need to be references. + alg_vars.append(instance.rr) + alg_vars.append(instance.L) + alg_vars.append(instance.V) + alg_vars.append(instance.FL) + + alg_eqns.append(instance.reflux_ratio) + alg_eqns.append(instance.flowrate_rectification) + alg_eqns.append(instance.vapor_column) + alg_eqns.append(instance.flowrate_stripping) + + input_vars = [pyo.Reference(instance.u1[:])] + + # Create a block to hold the reduced space model + reduced_space = pyo.Block(concrete=True) + reduced_space.obj = pyo.Reference(instance.REDUCED_SPACE_OBJ) + + n_input = len(input_vars) + + def differential_block_rule(b, i): + b.state = diff_vars[i] + b.deriv = deriv_vars[i] + b.disc = disc_eqns[i] + + def input_block_rule(b, i): + b.var = input_vars[i] + + reduced_space.differential_block = pyo.Block( + range(n_diff), + rule=differential_block_rule, + ) + reduced_space.input_block = pyo.Block( + range(n_input), + rule=input_block_rule, + ) + + reduced_space.external_block = ExternalGreyBoxBlock(instance.t) + + # Add reference to the constraint that specifies the initial conditions + reduced_space.init_rule = pyo.Reference(instance.init_rule) + + for t in instance.t: + if t == instance.t.first(): + reduced_space.external_block[t].deactivate() + continue + # Create and set external model for every external block + reduced_space_vars = ( + list(reduced_space.input_block[:].var[t]) + + list(reduced_space.differential_block[:].state[t]) + + list(reduced_space.differential_block[:].deriv[t]) + ) + external_vars = [v[t] for v in alg_vars] + residual_cons = [c[t] for c in diff_eqns] + external_cons = [c[t] for c in alg_eqns] + reduced_space.external_block[t].set_external_model( + ExternalPyomoModel( + reduced_space_vars, + external_vars, + residual_cons, + external_cons, + ), + inputs=reduced_space_vars, + ) + return reduced_space + + +def solve(reduced_space, tee=False): + solver = pyo.SolverFactory("cyipopt") + results = solver.solve(reduced_space, tee=tee) + return results diff --git a/performance/models/implicit/test_distill.py b/performance/models/implicit/test_distill.py new file mode 100644 index 0000000..cdf398f --- /dev/null +++ b/performance/models/implicit/test_distill.py @@ -0,0 +1,112 @@ +import os + +import pyomo.common.unittest as unittest + +from pyomo.common.timing import TicTocTimer +from pyomo.environ import * +from pyomo.dae import * + +_dir = os.path.dirname(__file__) + +from pyomo.contrib.pynumero.algorithms.solvers.cyipopt_solver import ( + cyipopt_available, +) + +@unittest.pytest.mark.performance +@unittest.skipUnless(cyipopt_available, "CyIpopt is not available") +class TestImplicitDistill(unittest.TestCase): + + _predicted_results = { + "u1[1]": 3.0, + "u1[2]": 1.0, + "u1[3]": 1.0, + "u1[4]": 1.0, + "u1[5]": 1.0, + "u1[6]": 1.0, + "u1[7]": 1.000000000894782, + "u1[8]": 1.0000001521950732, + "u1[9]": 1.1623645619362053, + "u1[10]": 1.3172203610288866, + "u1[11]": 1.4501249165622068, + "u1[12]": 1.5616919180444544, + "u1[13]": 1.6535269384632936, + "u1[14]": 1.7279135986010443, + "u1[15]": 1.7874071166551058, + "u1[16]": 1.8345254556710076, + "u1[17]": 1.8715664136378687, + "u1[18]": 1.900524442000351, + "u1[19]": 1.9230723526888294, + "u1[20]": 1.9405795753127466, + "u1[21]": 1.954147687752178, + "u1[22]": 1.9646514781090727, + "u1[23]": 1.9727790642832834, + "u1[24]": 1.9790679419304176, + "u1[25]": 1.9839358036166106, + "u1[26]": 1.987706049586084, + "u1[27]": 1.990628440391752, + "u1[28]": 1.9928955539673907, + "u1[29]": 1.9946557490664876, + "u1[30]": 1.9960232900793193, + "u1[31]": 1.9970862048604325, + "u1[32]": 1.997912354477815, + "u1[33]": 1.9985541054206861, + "u1[34]": 1.999051916788678, + "u1[35]": 1.9994370891585196, + "u1[36]": 1.9997338678771428, + "u1[37]": 1.9999610501879141, + "u1[38]": 2.000133211258838, + "u1[39]": 2.0002616372545323, + "u1[40]": 2.0003550326522737, + "u1[41]": 2.0004200528631024, + "u1[42]": 2.0004617009821075, + "u1[43]": 2.0004836185973205, + "u1[44]": 2.0004882949451854, + "u1[45]": 2.000477217021499, + "u1[46]": 2.0004509878661367, + "u1[47]": 2.000409457215106, + "u1[48]": 2.0003519538618404, + "u1[49]": 2.000277826301352, + "u1[50]": 2.0001878222716396, + "u1[51]": 2.000087833075182, + } + + def recordData(self, name, value): + """A method for recording data associated with a test. This method is only + meaningful when running this TestCase with 'nose', using the TestData plugin. + """ + tmp = getattr(self, 'testdata', None) + if not tmp is None: + tmp[name] = value + + def test_implicit_distill(self): + timer = TicTocTimer() + from .distill_DAE import ( + create_instance, + setup_implicit, + solve, + ) + instance = create_instance() + self.recordData('create_instance', timer.toc('create_instance')) + implicit_instance = setup_implicit(instance) + self.recordData('setup_implicit', timer.toc('setup_implicit')) + results = solve(implicit_instance) + self.recordData('solve', timer.toc('solve')) + + sol_data = { + name: instance.find_component(name).value + for name in self._predicted_results + } + self.assertStructuredAlmostEqual( + self._predicted_results, + sol_data, + delta=1e-5, + ) + self.recordData('check_results', timer.toc('check_results')) + + +if __name__ == '__main__': + import sys + from pyomo.common.fileutils import this_file_dir + sys.path.insert(0, os.path.dirname(this_file_dir())) + __package__ = os.path.basename(this_file_dir()) + unittest.main()