From bd1091f99ef0e6a7fa3f437b3f37bf973e425980 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Mon, 8 Aug 2022 17:39:26 -0400 Subject: [PATCH 1/6] init directory for implicit function --- performance/models/implicit/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 performance/models/implicit/__init__.py diff --git a/performance/models/implicit/__init__.py b/performance/models/implicit/__init__.py new file mode 100644 index 00000000..e69de29b From 8161a63c1f1315caa8b57a0765ea6e734736e8f5 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Mon, 8 Aug 2022 18:00:59 -0400 Subject: [PATCH 2/6] __doc__ for implicit.__init__ --- performance/models/implicit/__init__.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/performance/models/implicit/__init__.py b/performance/models/implicit/__init__.py index e69de29b..3bffb986 100644 --- a/performance/models/implicit/__init__.py +++ b/performance/models/implicit/__init__.py @@ -0,0 +1,21 @@ +__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. + + """ +) From dfcdbb3a28d78d3df6e892a5509b4198c60d88de Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Mon, 8 Aug 2022 18:15:43 -0400 Subject: [PATCH 3/6] update docstring --- performance/models/implicit/__init__.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/performance/models/implicit/__init__.py b/performance/models/implicit/__init__.py index 3bffb986..00f3bb9d 100644 --- a/performance/models/implicit/__init__.py +++ b/performance/models/implicit/__init__.py @@ -16,6 +16,16 @@ 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. """ ) From ae49a9a7f2d5d05c7c7ceb61946793013ff930a5 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Tue, 9 Aug 2022 14:56:23 -0400 Subject: [PATCH 4/6] initial files for distillation test problem --- performance/models/implicit/distill.dat | 12 ++ performance/models/implicit/distill_DAE.py | 159 ++++++++++++++++++++ performance/models/implicit/test_distill.py | 34 +++++ 3 files changed, 205 insertions(+) create mode 100644 performance/models/implicit/distill.dat create mode 100644 performance/models/implicit/distill_DAE.py create mode 100644 performance/models/implicit/test_distill.py diff --git a/performance/models/implicit/distill.dat b/performance/models/implicit/distill.dat new file mode 100644 index 00000000..c39307a8 --- /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 00000000..44dc2d1b --- /dev/null +++ b/performance/models/implicit/distill_DAE.py @@ -0,0 +1,159 @@ +# ___________________________________________________________________________ +# +# 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 diff --git a/performance/models/implicit/test_distill.py b/performance/models/implicit/test_distill.py new file mode 100644 index 00000000..46c8527b --- /dev/null +++ b/performance/models/implicit/test_distill.py @@ -0,0 +1,34 @@ +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__) + +@unittest.pytest.mark.performance +class TestStochPDEgas(unittest.TestCase): + + 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_stochpdegas_automatic(self): + timer = TicTocTimer() + from .distill_DAE import create_instance + instance = create_instance() + self.recordData('create_instance', timer.toc('create_instance')) + + +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() From 9e630ac4090c066f511d22e7692f208b5c9daeff Mon Sep 17 00:00:00 2001 From: Robert Date: Tue, 9 Aug 2022 16:02:48 -0400 Subject: [PATCH 5/6] implicit function DAE test --- performance/models/implicit/distill_DAE.py | 92 ++++++++++++++++++++- performance/models/implicit/test_distill.py | 79 +++++++++++++++++- 2 files changed, 167 insertions(+), 4 deletions(-) diff --git a/performance/models/implicit/distill_DAE.py b/performance/models/implicit/distill_DAE.py index 44dc2d1b..08cf576d 100644 --- a/performance/models/implicit/distill_DAE.py +++ b/performance/models/implicit/distill_DAE.py @@ -148,7 +148,8 @@ def rs_obj_rule(m): ) 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__) @@ -157,3 +158,92 @@ def create_instance(): 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 index 46c8527b..84791211 100644 --- a/performance/models/implicit/test_distill.py +++ b/performance/models/implicit/test_distill.py @@ -9,7 +9,61 @@ _dir = os.path.dirname(__file__) @unittest.pytest.mark.performance -class TestStochPDEgas(unittest.TestCase): +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 @@ -19,11 +73,30 @@ def recordData(self, name, value): if not tmp is None: tmp[name] = value - def test_stochpdegas_automatic(self): + def test_implicit_distill(self): timer = TicTocTimer() - from .distill_DAE import create_instance + 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__': From e05788c60c4f3ed841a1e14e0fdee4df8f047a64 Mon Sep 17 00:00:00 2001 From: Robert Date: Tue, 9 Aug 2022 16:32:11 -0400 Subject: [PATCH 6/6] skip test if cyipopt not available --- performance/models/implicit/test_distill.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/performance/models/implicit/test_distill.py b/performance/models/implicit/test_distill.py index 84791211..cdf398f8 100644 --- a/performance/models/implicit/test_distill.py +++ b/performance/models/implicit/test_distill.py @@ -8,7 +8,12 @@ _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 = {