diff --git a/doc/OnlineDocs/explanation/analysis/edi/additionaltips.rst b/doc/OnlineDocs/explanation/analysis/edi/additionaltips.rst new file mode 100644 index 00000000000..a09b2b830bd --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/edi/additionaltips.rst @@ -0,0 +1,37 @@ +Additional Tips +--------------- + +* Developers may need to install the following additional packages: + +:: + + pip install pytest + pip install pytest-cov + pip install sphinx + pip install sphinx_rtd_theme + pip install sphinx_copybutton + + +* If you wish to build the documentation locally, use: + +:: + + cd /pyomo/doc/OnlineDocs + make html + +then open the file ```` + + +* Unit tests and coverage can be run locally using: + +:: + + cd /pyomo/pyomo/contrib/edi + pytest --cov-report term-missing --cov=pyomo.contrib.edi -v ./tests/ + +or generating html output: + +:: + + cd /pyomo/pyomo/contrib/edi + pytest --cov-report html --cov=pyomo.contrib.edi -v ./tests/ diff --git a/doc/OnlineDocs/explanation/analysis/edi/advancedruntimeconstraints.rst b/doc/OnlineDocs/explanation/analysis/edi/advancedruntimeconstraints.rst new file mode 100644 index 00000000000..53a716c2f53 --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/edi/advancedruntimeconstraints.rst @@ -0,0 +1,115 @@ +Advanced Runtime Constraints +============================ + + +Automated unpacking for multi-use black-boxes +--------------------------------------------- + +Coding a black-box model often represents a significant effort, and it is therefore desirable to have the black-box work in many situations. A common case is to have a black-box model with a scalar input variable perform vectorized operations, ie, take in a vector and return a vector. In the more general case, this can be thought of as passing in multiple run-cases as input to the black-box model. + +The ``parseInputs()`` method enables this batch-like capability in a variety of forms: + +.. py:function:: BlackBoxFunctionModel.parseInputs(self, *args, **kwargs) + + Parses the inputs to a black-box model into a list of run-cases + + :param args: The self attribute in all python methods + :type self: black-box model + :param args: index passed arguments + :type args: list or tuple + :param kwargs: keyword passed arguments + :type kwargs: dict + + :return: runcases + :rtype: list + :return: returnMode + :rtype: int + :return: extras + :rtype: list + + +The function definition is not particularly helpful, so let's dive in a bit. For the typical user, we recommend that the top of all ``BlackBox()`` methods appear as follows: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 12 + :start-after: # BEGIN: AdvancedRTC_Snippet_02 + :end-before: # END: AdvancedRTC_Snippet_02 + + +Essentially, ``parseInputs()`` is a pre-processor that directly takes the inputs of the black-box. The ``parseInputs()`` method will check all of the inputs, ensure that size and units are correct, split into run cases as appropriate, and return a run-cases list that is ready to operate on. + +The ``runCases`` return (which can be named as any valid python name) is a list of dicts, where the keys to the dict are the names of the inputs declared in the ``__init__()`` method. Ex: ``runCases[0]['x']`` will give the ``'x'`` input (in units specified in the ``__init__()`` method) in the first run-case. + +The ``returnMode`` return (which can be named as any valid python name) is an integer that indicates how the return should be processed. If ``returnMode`` is passed as a ``kwarg``, then the passed in value will be cast to this output. If it is not passed in, then ``returnMode`` will be either ``-1*self.availableDerivative``, indicating that there is only a single run case, or ``self.availableDerivative`` which indicated multiple run cases. A negative value for ``returnMode`` indicates that there will be one less layer of indexing in the output. + +The ``extras`` return (which can be named as any valid python name) is a dict that will include all of the additional passed in values not expected by the black-box model. The extra ``kwargs`` will appear as normal, and extra ``args`` get put in a list in ``extras['remainingArgs']``. If you pass in a ``kwarg`` with key name ``'remainingArgs'``, it will be overwritten. The extras dict is the place to find options that may affect the code (ex: tolerances, run modes, etc) that are not expected inputs to the black-box model. + +Note that when using run case input, the output will always take the following unpacking structure: + +:: + + output[][0] = + output[][0][] = + + output[][1] = + output[][1][][] = + + +There is **not** a shortcut for single scalar output black boxes as is the case when not using run cases, the final index to output 0 must be included. + + + +An Example +++++++++++ + +There are many ways this functionality can be used, we provide an example here to get new users started + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: AdvancedRTC_Snippet_01 + :end-before: # END: AdvancedRTC_Snippet_01 + + +Check outputs +------------- + +There is a ``checkOutputs()`` method that is not supported in the current version. Contact the developers if you desire this functionality, but the following the practices described in this documentation should render the need for this moot. + + +Cases of non-scalar inputs or outputs +------------------------------------- + +Indexing can get somewhat complicated when inputs and outputs are not scalars. Users should be warned this feature is supported, but not well tested, so please contact the developers if you encounter any unusual behavior. + +The following unpacking remains the same: + +:: + + output[0] = + output[0][] = + + output[1] = + output[1][][] = + +However, for outputs, the result will be an array with dimensions equal to the size of the output. For Jacobians, it breaks down as the following: + +:: + + jacobian_list_entry[(output_dim_1_ix, output_dim_2_ix, ..., input_dim_1_ix, input_dim_2_ix, ...)] = + +For example, with an output that is 2x2 and an input that is also 2x2 + +:: + + output[1][][][(0,0,1,1)] + +is the derivative of ``output[0,0]`` with respect to ``input[1,1]`` + + + +Tips +---- + +* A model summary can be printed by calling ``print(model_instance.summary)`` diff --git a/doc/OnlineDocs/explanation/analysis/edi/blackboxconstraints.rst b/doc/OnlineDocs/explanation/analysis/edi/blackboxconstraints.rst new file mode 100644 index 00000000000..4d23f419354 --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/edi/blackboxconstraints.rst @@ -0,0 +1,411 @@ +Runtime (Black-Box) Constraints +=============================== + + +Overview +-------- + +While some constraints are explicitly known and can be written directly into the optimization problem, it is common (particularly in engineering design) for some relationships to be too complex to be directly coded as a constraint. + +EDI refers to these types of constraints as ``RuntimeConstraints`` because they are not constructed until they are needed by the solver. A particular subset of Runtime Constraints of interest are Black-Box constraints, that is, constraints which call to an external routine. To the average Pyomo and EDI user, ``RuntimeConstraints`` are (for all intents and purposes) Black-Box constraints, and the distinction is semantic. + +In other words, if you wish to code a black-box constraint using EDI, you will be using the Runtime Constraint constructor. + +In this context, a *Black-Box* is defined as a routine that performs hidden computation not visible EDI, Pyomo, or more generally the optimization algorithm. However, it is **not** assumed that black-boxes are unable to return gradient information. A black-box in this context may be capable of returning arbitrary derivative information. + + +Construction +------------ + +Runtime constraints consist of two separate elements that need to be constructed: an ``__init__`` function and a ``BlackBox`` function. Additionally, there are two pre-implemented functions that intermediate to advanced users will wish to interface with: ``BlackBox_Standardized`` and ``MultiCase``. + + +Constructing a Black Box +++++++++++++++++++++++++ + +First, we need to create an object which is visible to Pyomo/EDI that calls the black-box function. EDI calls this a ``BlackBoxFunctionModel``, and it is a base class that gets inherited into the objects you will create as a user. + +A simple example is shown below: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: RuntimeConstraints_Snippet_01 + :end-before: # END: RuntimeConstraints_Snippet_01 + +The inheriting classes can have any valid Python name (in this case ``Parabola``) and must have two methods ``__init__()`` and ``BlackBox()``. + + +The init method +*************** + +The ``__init__()`` function sets up the model, and has 5 distinct steps. First, the parent class ``__init__()`` must be called: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 16 + :start-after: # BEGIN: RuntimeConstraints_Snippet_02 + :end-before: # END: RuntimeConstraints_Snippet_02 + +In general, this line can be used verbatim. + +Next, you must tell the model what its inputs are by appending them to the ``self.inputs`` attribute. These inputs exist entirely in the local namespace of the black-box model, and are **independent** of the namespace in the optimization model (e.g. something called ``x`` in the optimization model can be called ``y`` in the black-box model). Inputs must have a ``name`` and ``units`` and can optionally have a ``description``, and ``size``. + +.. py:function:: self.inputs.append(name, units, description='', size=0) + + Appends a variable to a black box input list + + :param name: The name of the variable, any valid Python string. **Does not** have to match the name in the optimization formulation + :type name: str + :param units: The units of the variable. Every entry in a vector variable must have the same units. Entries of '', ' ', '-', 'None', and 'dimensionless' all become units.dimensionless. The units **must** be convertible to the units used in the optimization formulation (e.g., meters and feet), but **are not required** to be the same. Because of this, for example, a black-box can be written in imperial units, while an optimization formulation operates in metric. + :type units: str or pyomo.core.base.units_container._PyomoUnit + :param description: A description of the variable + :type description: str + :param size: The size (or shape) of the variable. Entries of 0, 1, and None all correspond to scalar variables. Other integers correspond to vector variables. Matrix and tensor variable are declared using lists of ints, ex: [10,10]. Matrix and tensor variables with a dimension of 1 (i.e., [10,10,1]) will be rejected as the extra dimension holds no meaningful value. + :type size: int or list + + +Models with multiple inputs simply call the ``self.input.append()`` command multiple times: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 16 + :start-after: # BEGIN: RuntimeConstraints_Snippet_03 + :end-before: # END: RuntimeConstraints_Snippet_03 + +Input names must be unique, and an error is raised if a repeated name is attempted to be set. + +Next, outputs must be added to the model. This is done identically to inputs, however the function is now ``self.outputs.append()`` + +.. py:function:: self.outputs.append(name, units, description='', size=0) + + Appends a variable to a black box output list + + :param name: The name of the variable, any valid Python string. **Does not** have to match the name in the optimization formulation + :type name: str + :param units: The units of the variable. Every entry in a vector variable must have the same units. Entries of '', ' ', '-', 'None', and 'dimensionless' all become units.dimensionless. The units **must** be convertible to the units used in the optimization formulation (ex, meters and feet), but **are not required** to be the same. Because of this, for example, a black-box can be written in imperial units, while an optimization formulation operates in metric. + :type units: str or pyomo.core.base.units_container._PyomoUnit + :param description: A description of the variable + :type description: str + :param size: The size (or shape) of the variable. Entries of 0, 1, and None all correspond to scalar variables. Other integers correspond to vector variables. Matrix and tensor variable are declared using lists of ints, ex: [10,10]. Matrix and tensor variables with a dimension of 1 (i.e., [10,10,1]) will be rejected as the extra dimension holds no meaningful value. + :type size: int or list + + +and similarly: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 16 + :start-after: # BEGIN: RuntimeConstraints_Snippet_04 + :end-before: # END: RuntimeConstraints_Snippet_04 + + +Finally, the highest available derivative must be set. For models being used in optimization, this will most often be ``1``, i.e. first derivative, gradient, or Jacobian information. + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 16 + :start-after: # BEGIN: RuntimeConstraints_Snippet_05 + :end-before: # END: RuntimeConstraints_Snippet_05 + + +The BlackBox method +******************* + +The ``BlackBox`` is extremely flexible, but here we present standard usage for a typical user. + +The ``BlackBox`` method should take in the inputs as arguments in the order defined during the ``__init__()`` method. Note that the method assumes inputs are passed in **with units** and outputs are returned **with units**. In general, the units on inputs and outputs need not be in any specific system, but should be convertible (e.g. meters and feet) to the units specified in the ``__init__()`` function. + +Various unpacking schemes are enabled by default via the ``parse_inputs`` function. Use of this function is not necessary, but provides for the parsing of index argumented lists (ex: ``function(x1, x2, x3)``) and keyword argumented dictionaries (ex: ``function({'x2':x2, 'x1':x1, 'x3',x3})``), along with a few other possibilities. + +The unit handling system in Pyomo can be rough at times, and so best practice is generally for the BlackBox function is expected to return values that are ``pint.Quantity`` types. These are obtained using the ``pyo.as_quantity()`` function. + +Since the units cannot be assumed on input, the first step in any black box is to convert to the model units: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 16 + :start-after: # BEGIN: RuntimeConstraints_Snippet_06 + :end-before: # END: RuntimeConstraints_Snippet_06 + + +And frequently, it is a good idea to cast these to a float value using ``pyomo.environ.value``: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 16 + :start-after: # BEGIN: RuntimeConstraints_Snippet_07 + :end-before: # END: RuntimeConstraints_Snippet_07 + +The assumed units can now be added if desired, but this may cause a slowdown in performance. Typical usage is to strip units then append at the end, unless many unit systems are being used in the actual computations. + +Operations can now be performed to compute the output and derivatives as desired. + +When preparing the outputs, note that all outputs must have units and be of type ``pint.Quantity``: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 16 + :start-after: # BEGIN: RuntimeConstraints_Snippet_08 + :end-before: # END: RuntimeConstraints_Snippet_08 + + +There are multiple options for packing the output. In general, the ``BlackBox`` method should output in a way that is convenient for the modeler. For simple black boxes, with less than 10 scalar inputs and scalar outputs, it is probably easiest to output as a tuple, as was done in the example here. Consider as a second example a function of form ``[u,v,w]=f(x,y,z)``. The simple packing would be: + +.. code-block:: python + + return [u,v,w], [[du_dx, du_dy, du_dz], [dv_dx, dv_dy, dv_dz], [dw_dx, dw_dy, dw_dz]] + +If only one output, the following is also allowed: + +.. code-block:: python + + return u, [du_dx, du_dy, du_dz] + +For more complex models, it is encouraged to switch to the ``NamedTuple`` output also used by the ``BlackBox_Standardized`` method: + +.. code-block:: python + + returnTuple = namedtuple('returnTuple', ['values', 'first', 'second']) + optTuple = namedtuple('optTuple', ['u','v','w']) + iptTuple = namedtuple('iptTuple', ['x','y','z']) + + values = optTuple(u,v,w) + first = optTuple( iptTuple(du_dx, du_dy, du_dz), + iptTuple(dv_dx, dv_dy, dv_dz), + iptTuple(dw_dx, dw_dy, dw_dz) ) + second = None # Second derivatives not currently supported + + return returnTuple(values,first,second) + +Dictionaries with the same keywords are also supported: + +.. code-block:: python + + values = {'u':u , 'v':v , 'w':w} + first = { 'u': {'x':du_dx ,'y':du_dy 'z':du_dz}, + 'v': {'x':dv_dx ,'y':dv_dy 'z':dv_dz}, + 'w': {'x':dw_dx ,'y':dw_dy 'z':dw_dz} } + second = None # Second derivatives not currently supported + + return { 'values': values, + 'first': first, + 'second': second } + +As are combinations of any of these options. + +In the event that the inputs and/or outputs are non-scalar, then outputs should be passed out as structures of their appropriate shape. Derivatives are a little more complicated. If the input **or** output is a scalar, then the derivative ``du_dx`` should have the shape of the non-scalar input/output. However, if **both** are non-scalar, then the output should be a numpy array that takes in indices of the inputs as indices. For example, an input (x) of dimension 2x2 and an output (u) of dimension 4x4x4, then the derivative information would be packed as: + +.. code-block:: python + + du_dx[0,0,0,0,0] # derivative of u[0,0,0] with respect to x[0,0] + du_dx[0,0,1,0,0] # derivative of u[0,0,1] with respect to x[0,0] + du_dx[0,0,0,1,1] # derivative of u[0,0,0] with respect to x[1,1] + +Note that this may change in the future, as developers are currently unsatisfied with extensions of this method to second order and higher derivatives. + +The BlackBox_Standardized Method +******************************** +The ``BlackBox`` method is designed to be highly flexible for use by practicing engineers. However, integrating these black box analysis tools into optimization often requires a common, structured framework for the output. The ``BlackBox_Standardized`` method provides this common interface. + +The ``BlackBox_Standardized`` method will **always** provide an output that is a nested series of ``NamedTuples``: + +- [0]--'values' + - [0]--'name_of_first_output': First output of the black box + - [1]--'name_of_first_output': Second output of the black box + - ... + - [n]--'name_of_last_output': Last output of the black box +- [1]--'first' + - [0]--'name_of_first_output' + - [0]--'name_of_first_input': Derivative of the first output of the black box wrt. the first input of the black box + - ... + - [n]--'name_of_last_input': Derivative of the first output of the black box wrt. the last input of the black box + - ... + - [n]--'name_of_last_output' + - [0]--'name_of_first_input': Derivative of the last output of the black box wrt. the first input of the black box + - ... + - [n]--'name_of_last_input': Derivative of the last output of the black box wrt. the last input of the black box +- [2]--'second' + - At present, this will always be ``None``, though second order support is planned in future versions + +For example, for a black box function ``[u,v,w] = f(x,y,z)``, the ``BlackBox_Standardized`` method would look as follows: + +.. code-block:: python + + opt = m.BlackBox_Standardized(x,y,z) + opt.values.u # The output u + opt[0][2] # The output w + opt.first.u.x # The derivative du/dx + opt.first.w.y # The derivative dw/dy + opt[1][1][2] # The derivative dv/dz + opt.second # Will be None in the current version + opt[2] # Will be None in the current version + +Note that while the current implementation does have a default ``BlackBox_Standardized`` method, this method is written for robustness and not performance. For users who want maximum performance, we recommend writing your own ``BlackBox_Standardized`` method and overriding the base class. Performance improvements will be most significant when a model has a large number of non-scalar inputs and/or outputs. + +Note also that the names in the named tuple are the names used in the black box model and **not** the names in the optimization problem. These two namespaces are intentionally separated. If you wish to use the namespace access, you must either get the names from ``[x.name for x in m.inputs]`` and ``[x.name for x in m.outputs]`` and use python's ``getattr`` function, or else have some other prior knowledge of the local black box namespace. + +Below is an example using the default ``BlackBox_Standardized`` method: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: RuntimeConstraints_Snippet_11 + :end-before: # END: RuntimeConstraints_Snippet_11 + +And now if we wish to implement our own custom ``BlackBox_Standardized`` method: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: RuntimeConstraints_Snippet_12 + :end-before: # END: RuntimeConstraints_Snippet_12 + +Note that if you are writing solvers for EDI, you should **exclusively** be calling the ``BlackBox_Standardized`` method, and never the basic ``BlackBox`` method, as this will ensure you are always working with a predictable data structure. + +As a modeler, you should **always** verify the output of ``BlackBox_Standardized`` as a final check in your effort to ensure it is packing the data as intended. Please notify the developers of any issues. + +The MultiCase Method +******************** +The ``MultiCase`` method provides a native capability to call the ``BlackBox`` method across multiple inputs simultaneously. This function is **not** vectorized in the base class and is **not** optimized for performance. If you wish to have a high performance vectorized function, you will need to implement your own method. + +Inputs to the ``MultiCase`` function should be a list of cases, which can be packed in any form accepted by the ``BlackBox_Standardized`` method. Overloading these functions may allow different forms of unpacking scheme. + +The output is a list of ``NamedTuple`` objects that are output from the ``BlackBox_Standardized`` method. If overloading, you may choose to output via a different packing scheme. + +Below is an example of overriding the default ``MultiCase`` method: + + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: RuntimeConstraints_Snippet_13 + :end-before: # END: RuntimeConstraints_Snippet_13 + + +Including a Black-Box in an EDI Formulation ++++++++++++++++++++++++++++++++++++++++++++ + +This second construction step is covered in the :doc:`Formulation <./formulation>` documentation, but is repeated here for completion. Future versions may differentiate this section. + +Runtime Constraints are declared one of two ways, just as regular constraints. The ``f.RuntimeConstraint()`` constructor is available: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_09 + :end-before: # END: Formulation_Snippet_09 + + +The ``f.RuntimeConstraint()`` constructor takes in the following inputs: + +.. py:function:: f.RuntimeConstraint(outputs, operators, inputs, black_box) + + Declares a runtime constraint in a pyomo.edi.formulation + + :param outputs: The outputs of the black box function + :type outputs: pyomo.environ.Var or list or tuple + :param operators: The operators that are used to construct constraints. Currently, only equality constraints are supported and will be the default no matter what is passed in here (see `this issue `__) + :type operators: str or list or tuple + :param inputs: The inputs to the black box function + :type inputs: pyomo.environ.Var or list or tuple + :param black_box: The object that stores the black-box function. See the :doc:`black box constraint documentation <./blackboxconstraints>` for details on constructing this object + :type black_box: pyomo.contrib.edi.BlackBoxFunctionModel + + +The following are alternative construction methods that may be of use: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_10 + :end-before: # END: Formulation_Snippet_10 + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_11 + :end-before: # END: Formulation_Snippet_11 + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_12 + :end-before: # END: Formulation_Snippet_12 + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_13 + :end-before: # END: Formulation_Snippet_13 + +However, more commonly we expect users to construct Runtime Constraints as a part of a ``f.ConstraintList()`` declaration. Simply include a list, tuple, or dict as a part of the ConstraintList as follows: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_14 + :end-before: # END: Formulation_Snippet_14 + +Any of the alternative declarations above are valid to pass into the ``f.ConstraintList()`` constructor, for example: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_15 + :end-before: # END: Formulation_Snippet_15 + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_16 + :end-before: # END: Formulation_Snippet_16 + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_17 + :end-before: # END: Formulation_Snippet_17 + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_18 + :end-before: # END: Formulation_Snippet_18 + + + +Examples +-------- + +More examples will be added over time. Feel free to reach out to the developers if you have questions regarding model development. + +A standard construction ++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: RuntimeConstraints_Snippet_10 + :end-before: # END: RuntimeConstraints_Snippet_10 + + +Tips +---- + +* Use the Pyomo ``tostr()`` function (``from pyomo.common.formatting import tostr``) to print the results of black-boxes for more meaningful printouts +* Align input and output declarations just as is recommended for optimization variable and constant declarations +* Declare an input/output all on one line, no matter what the style guides say +* * This interface is designed for subject matter experts who are not Python users to have a simple, easy path to include their tools/models into a Python based optimization architecture. +* Embrace units. They will save you so many times, it is well worth the minor additional overhead +* Pyomo units work slightly differently than pint (for those with pint experience), but those differences should be hidden from the model creator for the most part +* It is common to use this framework to call to a piece of software external to Python +* A model summary can be printed by calling ``print(model_instance.summary)`` + + +Known Issues +------------ + +* Currently only equality constraints are supported, pending an update to Pyomo (see `this issue `__) +* Runtime constraints must output to a variable, numbers and constants are not permitted (see `this issue `__) +* This functionality is not well tested when returning derivatives higher than first order. Though it should work, exercise caution and reach out to the dev team if questions arise. diff --git a/doc/OnlineDocs/explanation/analysis/edi/blackboxobjectives.rst b/doc/OnlineDocs/explanation/analysis/edi/blackboxobjectives.rst new file mode 100644 index 00000000000..cc2d869f72a --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/edi/blackboxobjectives.rst @@ -0,0 +1,6 @@ +Runtime (Black-Box) Objectives +============================== + +Black-box objectives are not currently supported and are awaiting an update to Pyomo (see `this issue `_). + +Please contact the developers if you are looking for this functionality diff --git a/doc/OnlineDocs/explanation/analysis/edi/constants.rst b/doc/OnlineDocs/explanation/analysis/edi/constants.rst new file mode 100644 index 00000000000..73e25234027 --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/edi/constants.rst @@ -0,0 +1,117 @@ +Constants +========= + +Overview +-------- +Constants are a key mechanism used to capture the relationship between variables. In engineering design, constants are often defined by physics or operational limits. + +The Constant constructor is a very thin wrapper on Pyomo :py:class:`Param `, and so experienced Pyomo users will not see any significant differences from base Pyomo. + + +Construction +------------ + +Constants are constructed by 1) creating an instance of a new constant in a EDI Formulation and 2) passing out this newly constructed constant to be used in objective and constraint construction. + +.. py:function:: f.Constant(name, value, units, description='', size=None, within=None) + + Declares a constant in a pyomo.contrib.edi.formulation + + :param name: The name of the constant for the purposes of tracking in the formulation. Commonly, this will be the same as the constant name in local namespace. + :type name: str + :param value: The value of the constant. For scalar constants, this should be a valid float or int for the specified domain. For vector constants, this will most often also be a single float or int, but a dictionary of index-value pairs is also accepted as in accordance with base Pyomo. NumPy arrays will be supported in a future release (see `this issue `__) + :type value: float or int or dict + :param units: The units of the constant. Every entry in a vector constant must have the same units. Entries of '', ' ', '-', 'None', and 'dimensionless' all become units.dimensionless + :type units: str or pyomo.core.base.units_container._PyomoUnit + :param description: A description of the constant + :type description: str + :param size: The size (or shape) of the constant. Entries of 0, 1, and None all correspond to scalar constants. Other integers correspond to vector constants. Matrix and tensor constants are declared using lists of ints, ex: [10,10]. Matrix and tensor constants with a dimension of 1 (i.e., [10,10,1]) will be rejected as the extra dimension holds no meaningful value. + :type size: int or list + :param within: The domain of the constant (ex: Reals, Integers, etc). Default of None constructs a constant in Reals. This option should rarely be used. + :type within: pyomo set + + :return: The constant that was declared in the formulation + :rtype: pyomo.core.base.param.ScalarParam or pyomo.core.base.param.IndexedParam + + +Relation to Pyomo Param +----------------------- + +The arguments ``name``, ``within``, and ``bounds`` are directly passed to the Pyomo ``Param`` constructor, with some minor checking. The ``value`` field is passed to the ``Param`` ``initialize`` field. The ``description`` field is passed to the ``doc`` field in the Pyomo ``Param``. Units are passed directly with an additional check. All Constants set the Pyomo ``Param`` ``mutable`` field to True. + +Non-scalar constants are constructed using Pyomo ``Sets``. Sets are constructed to be integer sets that fill the entire interval from lower bound to upper bound, ie a vector constant of length 5 would create a Pyomo ``Set`` with valid indices [0,1,2,3,4] with no skips. In this way, non-scalar constants are slightly less flexible than general non-scalar Pyomo ``Param``. + + +Examples +-------- + + +A standard declaration statement +++++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Constants_Snippet_01 + :end-before: # END: Constants_Snippet_01 + + +Shortest possible declaration ++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Constants_Snippet_02 + :end-before: # END: Constants_Snippet_02 + + +An alternative units definition ++++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Constants_Snippet_03 + :end-before: # END: Constants_Snippet_03 + + +A vector constant ++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Constants_Snippet_04 + :end-before: # END: Constants_Snippet_04 + + +A matrix/tensor constant +++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Constants_Snippet_05 + :end-before: # END: Constants_Snippet_05 + + +More complicated units definition ++++++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Constants_Snippet_06 + :end-before: # END: Constants_Snippet_06 + + +Tips +---- + +* Declare constants in alphabetical order. Trust me. It's a pain at first, but it saves a huge amount of time down the road, especially for large models. +* Designate a section in your file for constant declarations, as is done in the :doc:`introductory example <./quickstart>` +* Align all of your constant declarations in a pretty, grid like fashion. Depending on preference, these may or may not line up with variable declarations (I usually do not bother with this) +* Use the keyword names during constant declarations. Takes extra space, but is a massive boost to readability and intrepretability +* Declare one constant on one single line with no breaks, no matter what style guides tell you. Again, this is a significant boost to readability +* Do not skimp out on the description field, it is extremely helpful diff --git a/doc/OnlineDocs/explanation/analysis/edi/constraints.rst b/doc/OnlineDocs/explanation/analysis/edi/constraints.rst new file mode 100644 index 00000000000..bf283f93d31 --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/edi/constraints.rst @@ -0,0 +1,89 @@ +Constraints +=========== + + +Overview +-------- + +Constraints are the mathematical representation of rules that are imposed on your decisions/variables when minimizing or maximizing. In engineering design, constraints are often imposed by physics or operational limits. + +The Constraint constructor is a very thin wrapper on Pyomo :py:class:`Constraint `, and so experienced Pyomo users will not see any significant differences from base Pyomo. + + +Construction +------------ + +Constraints are constructed by creating an instance of a new constraint in a EDI Formulation + +.. py:function:: f.Constraint(expr) + + Declares a constraint in a pyomo.contrib.edi.formulation + + :param expr: The expression representing the constraint + :type expr: pyomo expression + + :return: None + :rtype: None + + +However, the expected use case is the ``f.ConstraintList()`` function: + +.. py:function:: f.ConstraintList(conList) + + Declares new constraints in a pyomo.contrib.edi.formulation from a list of inputs + + :param conList: The list of constraints to be generated. Entries will be Pyomo expressions, or lists/tuples/dicts that are used to create RuntimeConstraints (see :doc:`here <./blackboxconstraints>`) + :type conList: list + + :return: None + :rtype: None + + +Relation to Pyomo Constraint +---------------------------- + +The EDI constraint constructor is essentially a direct pass through to base Pyomo. Constraints will be added to the ``pyomo.ConcreteModel`` in increasing order with key ``constraint_###`` where the the index of the objective appears after the underscore. First constraint is labeled as ``constraint_1``, and constraint names are never padded with zeros. RuntimeConstraints also contribute to this counter. + + +Examples +-------- + +A standard declaration statement +++++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Constraints_Snippet_01 + :end-before: # END: Constraints_Snippet_01 + + +With the core constructor ++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Constraints_Snippet_02 + :end-before: # END: Constraints_Snippet_02 + +Using indexed variables and constants ++++++++++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Constraints_Snippet_03 + :end-before: # END: Constraints_Snippet_03 + +Tips +---- + +* For the typical user, constraints should always be constructed using the ``f.ConstraintList()`` function to produce a cleaner input file that is easier to modify + + +Known Issues +------------ + +* Indexed variables must be broken up using either indices or a Pyomo rule (see `this issue `__) +* Units that are inconsistent, but not the same (ie, meters and feet) will flag as invalid when checking units (see `this issue `__) diff --git a/doc/OnlineDocs/explanation/analysis/edi/examples.rst b/doc/OnlineDocs/explanation/analysis/edi/examples.rst new file mode 100644 index 00000000000..659f1b86de7 --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/edi/examples.rst @@ -0,0 +1,8 @@ +Examples +======== + +A Geometric Program for Aircraft Design +--------------------------------------- + +.. literalinclude:: ../../../../pyomo/contrib/edi/examples/aircraft_gp.py + :language: python diff --git a/doc/OnlineDocs/explanation/analysis/edi/formulation.rst b/doc/OnlineDocs/explanation/analysis/edi/formulation.rst new file mode 100644 index 00000000000..64968eff425 --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/edi/formulation.rst @@ -0,0 +1,206 @@ +Formulation +=========== + +.. |br| raw:: html + +
+ + +The core object in EDI is called a *Formulation*. The core object in EDI is called a *Formulation*. For experienced Pyomo users, a Formulation inherits from a Pyomo :class:`~pyomo.environ.ConcreteModel`, and can therefore be treated exactly as a typical :class:`~pyomo.environ.ConcreteModel` with a few additional features. + +Each modeling element (ex: Variable, Constant, Objective, and Constraint) has a constructor that is used to create the corresponding element in the *Formulation* instance. In addition, there are a number of helper functions that collect and return model elements or perform supporting actions. + +Construction +------------ +A *Formulation* is constructed as follows: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_01 + :end-before: # END: Formulation_Snippet_01 + + +Standard practice is to construct a formulation to namespace variable ``f``, but any valid Python name can be used. Standard Pyomo practice would be to construct this to ``model`` or ``m``. + + +Declaring Variables +------------------- + +See the :doc:`Variables <./variables>` Documentation + +Variables are declared using the ``f.Variable()`` function. This function creates an instance of :class:`~pyomo.environ.Var` and adds it to the ``pyomo.contrib.edi.Formulation``. The function returns an instance of a :class:`~pyomo.environ.Var` that can be used in later construction. + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_02 + :end-before: # END: Formulation_Snippet_02 + + +Declaring Constants +------------------- + +See the :doc:`Constants <./constants>` Documentation + +Constants (referred to in base Pyomo as parameters or ``Params``) are declared using the ``f.Constant()`` function. This function creates an instance of a :class:`~pyomo.environ.Param` and adds it to the ``pyomo.contrib.ediFormulation``. This function also returns an instance of :class:`~pyomo.environ.Param` that can be used in later construction. + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_03 + :end-before: # END: Formulation_Snippet_03 + + + +Declaring Objectives +-------------------- + +See the :doc:`Objectives <./objectives>` Documentation + +Objectives are declared using the ``f.Objective()`` function. This function creates an instance of a :class:`~pyomo.environ.Objective` and adds it to the ``pyomo.contrib.edi.Formulation``. Multiple objectives can be declared, but interpretation of multiple objectives will depend on the solver. The returned values of the ``f.Variable()`` and ``f.Constant()`` declarations can be used to construct the objective. Black-box (i.e. Runtime) objectives are not supported at this time, but are planned in a future update. + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_04 + :end-before: # END: Formulation_Snippet_04 + + +By default, objectives are minimized, but can be switched to a maximize using the ``sense`` keyword from Pyomo: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_05 + :end-before: # END: Formulation_Snippet_05 + + +Note: Future version will allow a string to be passed into ``sense`` (see `this issue `_) + + +Declaring Constraints +--------------------- + +See the :doc:`Constraints <./constraints>` Documentation + +Constraints can be declared in two ways. First is using the standard ``f.Constraint()`` constructor. This function creates an instance of :class:`~pyomo.environ.Constraint` and adds it to the ``pyomo.contrib.edi.Formulation``. The operators ``<=``, ``>=``, and ``==`` are used as constraint constructors. + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_06 + :end-before: # END: Formulation_Snippet_06 + +Constraints can also be declared using the ``f.ConstraintList()`` function. This function takes in a list of constraints and allows for multiple constraints to be declared in one go: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_07 + :end-before: # END: Formulation_Snippet_07 + +The constraint list can also be declared a priori and passed in to the ``f.ConstraintList()`` function, which may be beneficial for complex models: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_08 + :end-before: # END: Formulation_Snippet_08 + +We recommend that most users should be using the ``f.ConstraintList()`` function, with the ``f.Constraint()`` function being reserved for under-the-hood usage. + + +Declaring Runtime (Black-Box) Constraints +----------------------------------------- + +See the :doc:`Runtime (Black-Box) Constraints <./blackboxconstraints>` Documentation + +One of the main features of EDI is the streamlined implementation of Black-Box Constraints. A *Black-Box* is defined as a routine that performs hidden computation not visible to EDI, Pyomo, or more generally the optimization algorithm. However, it is **not** assumed that black-boxes are unable to return gradient information. A black-box in this context may be capable of returning arbitrary derivative information. + +Black-box constraints are considered to be a sub-class of a more general class of constraints called *Runtime Constraints*, that is constraints that are not actually constructed until the optimization routine is running. In most cases, Runtime Constraints are approximated as linear by the solver, and therefore a Runtime Constraint is expected to provide function evaluations and gradient information. + +The use of Runtime Constraints requires a black box model that is discussed in detail in the dedicated documentation (see :doc:`here <./blackboxconstraints>`), but for the purposes of demonstrating the constructors, a simple black box example will appear in all of the code snippets below. + +Runtime Constraints are declared one of two ways, just as regular constraints. The ``f.RuntimeConstraint()`` constructor is available: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_09 + :end-before: # END: Formulation_Snippet_09 + + +The ``f.RuntimeConstraint()`` constructor takes in the following inputs: + +.. py:function:: f.RuntimeConstraint(outputs, operators, inputs, black_box) + :noindex: + + Declares a runtime constraint in a pyomo.edi.formulation + + :param outputs: The outputs of the black box function + :type outputs: pyomo.environ.Var or list or tuple + :param operators: The operators that are used to construct constraints. Currently, only equality constraints are supported and will be the default no matter what is passed in here (see `this issue `__) + :type operators: str or list or tuple + :param inputs: The inputs to the black box function + :type inputs: pyomo.environ.Var or list or tuple + :param black_box: The object that stores the black-box function. See the :doc:`black box constraint documentation <./blackboxconstraints>` for details on constructing this object + :type black_box: pyomo.contrib.edi.BlackBoxFunctionModel + +However, more commonly we expect users to construct Runtime Constraints as a part of a ``f.ConstraintList()`` declaration. Simply include a list, tuple, or dict as a part of the ConstraintList as follows: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_14 + :end-before: # END: Formulation_Snippet_14 + +The ``f.ConstraintList()`` constructor tests for types ``tuple``, ``list``, or ``dict``, and attempts to construct a ``RuntimeConstraint`` from this data. Types of ``tuple`` and ``list`` should include 4 elements that are the ordered entries to the ``RuntimeConstraint`` constructor (ie, ``*args``), and type ``dict`` should have keyword arguments to the ``RuntimeConstraint`` constructor (ie, ``**kwrgs``). Examples are as follows: + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_15 + :end-before: # END: Formulation_Snippet_15 + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_16 + :end-before: # END: Formulation_Snippet_16 + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_17 + :end-before: # END: Formulation_Snippet_17 + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Formulation_Snippet_18 + :end-before: # END: Formulation_Snippet_18 + +Support Functions +----------------- + +``f.get_variables()`` |br| +Returns a list variables that have been defined in the formulation in declaration order. Will only return variables defined via EDI. + +``f.get_constants()`` |br| +Returns a list of constants that have been defined in the formulation in declaration order. Will only return constants/parameters defined via EDI. + +``f.get_objectives()`` |br| +Returns a list of objectives that have been defined in the formulation in declaration order. Will only return objectives defined via EDI. + +``f.get_constraints()`` |br| +Returns a list of constraints that have been defined in the formulation in declaration order. This command returns a list that includes both explicit and runtime (black-box) constraints, but only constraints that have been defined via EDI. + +``f.get_explicitConstraints()`` |br| +Returns a list of *explicit* constraints that have been defined in the formulation in declaration order. This command returns a list that includes *only* the explicit constraints and *not* the runtime (black-box) constraints. Only includes constraints that have been defined via EDI. + +``f.get_runtimeConstraints()`` |br| +Returns a list of *runtime* (ie. black-box) constraints that have been defined in the formulation in declaration order. This command returns a list that includes *only* the runtime constraints and *not* the explicit constraints. Only includes constraints that have been defined via EDI. + +``f.check_units()`` |br| +Checks the units of each objective and constraint for consistency. Will only check objectives and constraints defined via EDI. diff --git a/doc/OnlineDocs/explanation/analysis/edi/index.rst b/doc/OnlineDocs/explanation/analysis/edi/index.rst new file mode 100644 index 00000000000..244262a64c4 --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/edi/index.rst @@ -0,0 +1,47 @@ +Engineering Design Interface (EDI) +================================== + +The Pyomo Engineering Design Interface (EDI) is a lightweight wrapper on the Pyomo language that is targeted at composing engineering design optimization problems. The language and interface have been designed to mimic many of the features found in `GPkit `_ and `CVXPY `_ while also providing a simple, clean interface for black-box analysis codes that are common in engineering design applications. + + +Installation +------------ + +EDI installs as a part of the standard Pyomo install: + +:: + + pip install pyomo + + +EDI also requires packages that are optional in base Pyomo: + +:: + + pip install numpy + pip install scipy + pip install pint + + +User's Guide +------------ + +.. toctree:: + :maxdepth: 4 + + quickstart.rst + formulation.rst + variables.rst + constants.rst + objectives.rst + blackboxobjectives.rst + constraints.rst + blackboxconstraints.rst + examples.rst + additionaltips.rst + + +Developers +---------- + +The Pyomo EDI interface is developed and maintained by `Cody Karcher `_ diff --git a/doc/OnlineDocs/explanation/analysis/edi/objectives.rst b/doc/OnlineDocs/explanation/analysis/edi/objectives.rst new file mode 100644 index 00000000000..9144fbd42e5 --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/edi/objectives.rst @@ -0,0 +1,107 @@ +Objectives +========== + +Overview +-------- + +Objectives are the mathematical representation of what you wish to minimize or maximize. In engineering design, objectives are often minimizing cost, material, or time, or alternatively maximizing profit or utility. + +The Objective constructor is a very thin wrapper on Pyomo :py:class:`Objective `, and so experienced Pyomo users will not see any significant differences from base Pyomo. + + +Construction +------------ + +Objectives are constructed by creating an instance of a new objective in a EDI Formulation + +.. py:function:: f.Objective(expr, sense=minimize) + + Declares an objective in a pyomo.contrib.edi.formulation + + :param expr: The expression to be optimized + :type expr: pyomo expression + :param sense: The sense in which the objective should be optimized, either minimized or maximized. Can import ``minimize`` and ``maximize`` from ``pyomo.environ``, but minimize corresponds to an integer of 1 and maximize to an integer of -1. + :type sense: int + + :return: None + :rtype: None + + +Relation to Pyomo Objective +--------------------------- + +The EDI objective constructor is essentially a direct pass through to base Pyomo. Objectives will be added to the ``pyomo.ConcreteModel`` in increasing order with key ``objective_###`` where the the index of the objective appears after the underscore. First objective is labeled as ``objective_1``, and objective names are never padded with zeros. + + +Examples +-------- + +A standard declaration statement +++++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Objectives_Snippet_01 + :end-before: # END: Objectives_Snippet_01 + +With a non-linear objective ++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Objectives_Snippet_02 + :end-before: # END: Objectives_Snippet_02 + +Explicitly minimize ++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Objectives_Snippet_03 + :end-before: # END: Objectives_Snippet_03 + +Explicitly minimize using integer ++++++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Objectives_Snippet_04 + :end-before: # END: Objectives_Snippet_04 + + +Maximizing +++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Objectives_Snippet_05 + :end-before: # END: Objectives_Snippet_05 + +Maximizing using integer +++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Objectives_Snippet_06 + :end-before: # END: Objectives_Snippet_06 + +Using indexed variables and constants ++++++++++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Objectives_Snippet_07 + :end-before: # END: Objectives_Snippet_07 + + +Tips +---- + +* Objectives are a pretty natural place to break your file. Put at least one blank line above and below the objective constructor and use good sectioning to create a whitespace easily identifiable when scrolling quickly diff --git a/doc/OnlineDocs/explanation/analysis/edi/quickstart.rst b/doc/OnlineDocs/explanation/analysis/edi/quickstart.rst new file mode 100644 index 00000000000..cb90b4e93ef --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/edi/quickstart.rst @@ -0,0 +1,22 @@ +Quickstart +========== + +Below is a simple example written using EDI that should help get new users off the ground. For simple problems this example case is fine to build on, but for more advanced techniques (ex: using multi-dimensional variables or advanced black-box usage) see the additional documentation. + +.. note:: + + The files and examples that appear in this documentation have been adapted to conform to PEP-8 standards, however, we **strongly discourage** the adoption of PEP-8 in the course of regular EDI usage as it makes the code (particularly variable declarations, constant declarations, and constraint list declarations) extremely challenging to read and debug + +The example shown here minimizes a linear objective function subject to the interior area of the unit circle: + +.. math:: + \begin{align*} + & \underset{x,y,z}{\text{minimize}} + & & x+y \\ + & \text{subject to} + & & z = x^2 + y^2\\ + &&& z \leq 1.0 \text{ [m$^2$]} + \end{align*} + +.. literalinclude:: ../../../../pyomo/contrib/edi/examples/readme_example.py + :language: python diff --git a/doc/OnlineDocs/explanation/analysis/edi/variables.rst b/doc/OnlineDocs/explanation/analysis/edi/variables.rst new file mode 100644 index 00000000000..87c584852eb --- /dev/null +++ b/doc/OnlineDocs/explanation/analysis/edi/variables.rst @@ -0,0 +1,140 @@ +Variables +========= + +Overview +-------- +Variables are the mathematical representation of individual decisions being considered by the optimizer. In engineering design, variables are often representations of geometric parameters and operating conditions. + +Using the EDI package, variables can be defined as both scalar (``pyomo.core.base.var.ScalarVar``) and vector/matrix/tensor (``pyomo.core.base.var.IndexedVar``), and can exist in many mathematical spaces (All Real, Integers, etc). + +The Variable constructor is a very thin wrapper, and so experienced Pyomo users will not see any significant differences from base Pyomo. + + +Construction +------------ + +Variables are constructed by 1) creating an instance of a new variable in a EDI Formulation and 2) passing out this newly constructed variable to be used in objective and constraint construction. + +.. py:function:: f.Variable(name, guess, units, description='', size=None, bounds=None, domain=None) + + Declares a variable in a pyomo.contrib.edi.formulation + + :param name: The name of the variable for the purposes of tracking in the formulation. Commonly, this will be the same as the variable name in local namespace. + :type name: str + :param guess: The initial guess of the variable. For scalar variables, this should be a valid float or int for the specified domain. For vector variables, this will most often also be a single float or int, but a dictionary of index-value pairs is also accepted as in accordance with base Pyomo. NumPy arrays will be supported in a future release (see `this issue `_) + :type guess: float or int or dict + :param units: The units of the variable. Every entry in a vector variable must have the same units. Entries of '', ' ', '-', 'None', and 'dimensionless' all become units.dimensionless + :type units: str or pyomo.core.base.units_container._PyomoUnit + :param description: A description of the variable + :type description: str + :param size: The size (or shape) of the variable. Entries of 0, 1, and None all correspond to scalar variables. Other integers correspond to vector variables. Matrix and tensor variable are declared using lists of ints, ex: [10,10]. Matrix and tensor variables with a dimension of 1 (ie, [10,10,1]) will be rejected as the extra dimension holds no meaningful value. + :type size: int or list + :param bounds: The bounds on the variable. A list or tuple of two elements [lower_bound, upper_bound] where the two bounds are assumed to be either ints or floats. WARNING: User is currently responsible for ensuring the units are correct (see `this issue `__) + :type bounds: list or tuple + :param domain: The domain of the variable (ex: Reals, Integers, etc). Default of None constructs a variable in Reals. + :type domain: pyomo set + + :return: The variable that was declared in the formulation + :rtype: pyomo.core.base.var.ScalarVar or pyomo.core.base.var.IndexedVar + + +Relation to Pyomo Var +--------------------- + +The fields: name, domain, and bounds are directly passed to the Pyomo ``Var`` constructor, with some minor checking. The guess field is passed to initialize. The description field is passed to the doc field in the Pyomo ``Var``. Units are passed directly with an additional check. + +Non-scalar variables are constructed using Pyomo ``Sets``. Sets are constructed to be integer sets that fill the entire interval from lower bound to upper bound, ie a vector variable of length 5 would create a Pyomo ``Set`` with valid indices [0,1,2,3,4] with no skips. In this way, non-scalar constants are slightly less flexible than general non-scalar Pyomo ``Params``. + + +Examples +-------- + + +A standard declaration statement +++++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Variables_Snippet_01 + :end-before: # END: Variables_Snippet_01 + + +Shortest possible declaration ++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Variables_Snippet_02 + :end-before: # END: Variables_Snippet_02 + + +A variable with bounds +++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Variables_Snippet_03 + :end-before: # END: Variables_Snippet_03 + + +An integer variable ++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Variables_Snippet_04 + :end-before: # END: Variables_Snippet_04 + +An alternative units definition ++++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Variables_Snippet_05 + :end-before: # END: Variables_Snippet_05 + + +A vector variable ++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Variables_Snippet_06 + :end-before: # END: Variables_Snippet_06 + + +A matrix/tensor variable +++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Variables_Snippet_07 + :end-before: # END: Variables_Snippet_07 + + +More complicated units definition ++++++++++++++++++++++++++++++++++ + +.. literalinclude:: ../../../../pyomo/contrib/edi/tests/test_docSnippets.py + :language: python + :dedent: 8 + :start-after: # BEGIN: Variables_Snippet_08 + :end-before: # END: Variables_Snippet_08 + + +Tips +---- + +* Declare variables in alphabetical order. Trust me. It's a pain at first, but it saves a huge amount of time down the road, especially for large models. +* Designate a section in your file for variable declarations, as is done in the :doc:`introductory example <./quickstart>` +* Align all of your variable declarations in a pretty, grid like fashion. Again, a pain at first, but a big time saver +* Use the keyword names during variable declarations. Takes extra space, but is a massive boost to readability and intrepretability +* Declare one variable on one single line with no breaks, no matter what style guides tell you. Again, this is a significant boost to readability +* Do not skimp out on the description field, it is extremely helpful diff --git a/doc/OnlineDocs/explanation/analysis/index.rst b/doc/OnlineDocs/explanation/analysis/index.rst index 0a8e3c3b416..a7899b29b4b 100644 --- a/doc/OnlineDocs/explanation/analysis/index.rst +++ b/doc/OnlineDocs/explanation/analysis/index.rst @@ -7,6 +7,7 @@ Analysis in Pyomo alternative_solutions community doe/doe + edi/index iis incidence/index mpc/index diff --git a/doc/OnlineDocs/explanation/contrib_index.txt b/doc/OnlineDocs/explanation/contrib_index.txt deleted file mode 100644 index 65c14a721df..00000000000 --- a/doc/OnlineDocs/explanation/contrib_index.txt +++ /dev/null @@ -1,47 +0,0 @@ -Third-Party Contributions -========================= - -Pyomo includes a variety of additional features and functionality -provided by third parties through the ``pyomo.contrib`` package. This -package includes both contributions included with the main Pyomo -distribution and wrappers for third-party packages that must be -installed separately. - -These packages are maintained by the original contributors and are -managed as *optional* Pyomo packages. - -Contributed packages distributed with Pyomo: - -.. toctree:: - :maxdepth: 1 - - alternative_solutions.rst - community.rst - doe/doe.rst - gdpopt.rst - iis.rst - incidence/index.rst - latex_printer.rst - mindtpy.rst - mpc/index.rst - multistart.rst - preprocessing.rst - parmest/index.rst - pynumero/index.rst - pyros.rst - sensitivity_toolbox.rst - trustregion.rst - -Contributed Pyomo interfaces to other packages: - -.. toctree:: - :maxdepth: 1 - - mcpp.rst - satsolver.rst - - -Contributed packages distributed independently of Pyomo, but accessible -through ``pyomo.contrib``: - -* `pyomo.contrib.simplemodel `_ diff --git a/pyomo/contrib/edi/README.md b/pyomo/contrib/edi/README.md new file mode 100644 index 00000000000..7982834f4f6 --- /dev/null +++ b/pyomo/contrib/edi/README.md @@ -0,0 +1,101 @@ +# Engineering Design Interface + +The Pyomo Engineering Design Interface (EDI) is a lightweight wrapper on the Pyomo language that is targeted at composing engineering design optimization problems. The language and interface have been designed to mimic many of the features found in [GPkit](https://github.com/convexengineering/gpkit) and [CVXPY](https://github.com/cvxpy/cvxpy) while also providing a simple, clean interface for black-box analysis codes that are common in engineering design applications. + +## Installation + +EDI is a part of the standard installation process for Pyomo: +``` +pip install pyomo +``` + +EDI also requires the pint dependency that is optional in base Pyomo: +``` +pip install pint +``` + +## Usage + +The core object in EDI is the `Formulation` object, which inherits from the `pyomo.environ.ConcreteModel`. Essentially, a `Formulation` is a Pyomo `Model` with some extra stuff, but can be treated exactly as if it were a Pyomo `Model`. However, an EDI `Formulation` has some additional features that can help simplify model construction. + +```python +# ================= +# Import Statements +# ================= +import pyomo.environ as pyo +from pyomo.environ import units +from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + +# =================== +# Declare Formulation +# =================== +f = Formulation() + +# ================= +# Declare Variables +# ================= +x = f.Variable(name = 'x', guess = 1.0, units = 'm' , description = 'The x variable') +y = f.Variable(name = 'y', guess = 1.0, units = 'm' , description = 'The y variable') +z = f.Variable(name = 'z', guess = 1.0, units = 'm^2', description = 'The unit circle output') + +# ================= +# Declare Constants +# ================= +c = f.Constant(name = 'c', value = 1.0, units = '', description = 'A constant c', size = 2) + +# ===================== +# Declare the Objective +# ===================== +f.Objective( + c[0]*x + c[1]*y +) + +# =================== +# Declare a Black Box +# =================== +class UnitCircle(BlackBoxFunctionModel): + def __init__(self): # The initialization function + + # Initialize the black box model + super().__init__() + + # A brief description of the model + self.description = 'This model evaluates the function: z = x**2 + y**2' + + # Declare the black box model inputs + self.inputs.append(name = 'x', units = 'ft' , description = 'The x variable') + self.inputs.append(name = 'y', units = 'ft' , description = 'The y variable') + + # Declare the black box model outputs + self.outputs.append(name = 'z', units = 'ft**2', description = 'Resultant of the unit circle evaluation') + + # Declare the maximum available derivative + self.availableDerivative = 1 + + # Post-initialization setup + self.post_init_setup() + + def BlackBox(self, x, y): # The actual function that does things + x = pyo.value(units.convert(x,self.inputs['x'].units)) # Converts to correct units then casts to float + y = pyo.value(units.convert(y,self.inputs['y'].units)) # Converts to correct units then casts to float + + z = x**2 + y**2 # Compute z + dzdx = 2*x # Compute dz/dx + dzdy = 2*y # Compute dz/dy + + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + +# ======================= +# Declare the Constraints +# ======================= +f.ConstraintList( + [ + [ z, '==', [x,y], UnitCircle() ] , + z <= 1*units.m**2 + ] +) +``` diff --git a/pyomo/contrib/edi/__init__.py b/pyomo/contrib/edi/__init__.py new file mode 100644 index 00000000000..6c38af6220c --- /dev/null +++ b/pyomo/contrib/edi/__init__.py @@ -0,0 +1,80 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# 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. +# +# Development of this module was conducted as part of the Institute for +# the Design of Advanced Energy Systems (IDAES) with support through the +# Simulation-Based Engineering, Crosscutting Research Program within the +# U.S. Department of Energy’s Office of Fossil Energy and Carbon Management. +# +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +# # Register the pint quantity type to prevent warnings +# from pyomo.common.numeric_types import RegisterNumericType +# # try: +# import pint +# RegisterNumericType(pint.Quantity) +# # except: +# # pass + +# Recommended just to build all of the appropriate things +import pyomo.environ as pyo + +# Import the relevant classes from Formulation +try: + from pyomo.contrib.edi.formulation import Formulation +except: + pass + # in this case, the dependencies are not installed, nothing will work + +# Import the black box modeling tools +try: + from pyomo.contrib.edi.blackBoxFunctionModel import BlackBoxFunctionModel + +except: + pass + # in this case, the dependencies are not installed, nothing will work + + +# the printer that does units ok +import copy +from pyomo.core.base.units_container import _PyomoUnit +from pyomo.core.expr.numeric_expr import NPV_ProductExpression, NPV_DivisionExpression +from collections import namedtuple + +from pyomo.common.dependencies import numpy, numpy_available + +if numpy_available: + np = numpy + + def recursive_sub(x_in): + x = list(copy.deepcopy(x_in)) + for i in range(0, len(x)): + if isinstance(x[i], _PyomoUnit): + x[i] = '1.0*' + str(x[i]) + elif ( + isinstance( + x[i], (NPV_ProductExpression, NPV_DivisionExpression, np.float64) + ) + or x[i] is None + ): + if pyo.value(x[i]) == 1: + x[i] = '1.0*' + str(x[i]) + else: + x[i] = str(x[i]) + else: + x[i] = recursive_sub(list(x[i])) + return x + + def ediprint(x): + print(recursive_sub(x)) + +else: + pass + # in this case, the dependencies are not installed, nothing will work diff --git a/pyomo/contrib/edi/blackBoxFunctionModel.py b/pyomo/contrib/edi/blackBoxFunctionModel.py new file mode 100644 index 00000000000..b8677bf406f --- /dev/null +++ b/pyomo/contrib/edi/blackBoxFunctionModel.py @@ -0,0 +1,1409 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# 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. +# +# Development of this module was conducted as part of the Institute for +# the Design of Advanced Energy Systems (IDAES) with support through the +# Simulation-Based Engineering, Crosscutting Research Program within the +# U.S. Department of Energy’s Office of Fossil Energy and Carbon Management. +# +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import copy +import itertools +from collections import namedtuple +import pyomo +import pyomo.environ as pyo +from pyomo.environ import units as pyomo_units +from pyomo.core.base.units_container import as_quantity +from pyomo.common.dependencies import attempt_import + +from pyomo.core.expr.ndarray import NumericNDArray + +# from pyomo.common.numeric_types import RegisterNumericType +try: + import pint +except: + raise ImportError( + "pyomo.contrib.edi requires pint to enable black box capability, fix with 'pip install pint' " + ) + + +egb, egb_available = attempt_import( + "pyomo.contrib.pynumero.interfaces.external_grey_box" +) + +from pyomo.common.dependencies import numpy, numpy_available +from pyomo.common.dependencies import scipy, scipy_available + + +if numpy_available: + np = numpy +else: + raise ImportError( + "pyomo.contrib.edi requires numpy to enable black box capability, fix with 'pip install numpy' " + ) + +if scipy_available: + sps = scipy.sparse +else: + raise ImportError( + "pyomo.contrib.edi requires scipy to enable black box capability, fix with 'pip install scipy' " + ) + +if egb_available: + from pyomo.contrib.pynumero.interfaces.external_grey_box import ( + ExternalGreyBoxModel, + ExternalGreyBoxBlock, + ) +else: + raise ImportError( + "pyomo.contrib.edi requires pyomo.contrib.pynumero to be installed to enable black box capability, this should have installed with base pyomo" + ) + + +class BlackBoxFunctionModel_Variable(object): + def __init__(self, name, units, description='', size=0): + # Order matters + self.name = name + self.units = units + self.size = size + self.description = description + + # ===================================================================================================================== + # The printing function + # ===================================================================================================================== + def __repr__(self): + return self.name + + # ===================================================================================================================== + # Define the name + # ===================================================================================================================== + @property + def name(self): + return self._name + + @name.setter + def name(self, val): + if isinstance(val, str): + self._name = val + else: + raise ValueError('Invalid name. Must be a string.') + + # ===================================================================================================================== + # Define the units + # ===================================================================================================================== + @property + def units(self): + return self._units + + @units.setter + def units(self, val): + # set dimensionless if a null string is passed in + if isinstance(val, str): + if val in ['-', '', ' ']: + val = 'dimensionless' + if val is None: + val = 'dimensionless' + + if isinstance(val, str): + self._units = pyomo_units.__getattr__(val) + elif isinstance(val, pyomo.core.base.units_container._PyomoUnit): + self._units = val + else: + raise ValueError( + 'Invalid units. Must be a string compatible with pint or a unit instance.' + ) + + # ===================================================================================================================== + # Define the size + # ===================================================================================================================== + @property + def size(self): + return self._size + + @size.setter + def size(self, val): + invalid = False + if isinstance(val, (list, tuple)): + sizeTemp = [] + for x in val: + if isinstance(x, str): + # is a vector of unknown length, should be 'inf', but any string accepted + x = -1 + # pass + elif not isinstance(x, int): + raise ValueError( + 'Invalid size. Must be an integer or list/tuple of integers' + ) + if x == 1: + raise ValueError( + 'A value of 1 is not valid for defining size. Use fewer dimensions.' + ) + sizeTemp.append(x) + self._size = sizeTemp + else: + if val is None: + self._size = 0 # set to scalar + elif isinstance(val, str): + # is a 1D vector of unknown length, should be 'inf', but any string accepted + self._size = -1 + # pass + elif isinstance(val, int): + if val == 1: + raise ValueError( + 'A value of 1 is not valid for defining size. Use 0 to indicate a scalar value.' + ) + else: + self._size = val + else: + raise ValueError( + 'Invalid size. Must be an integer or list/tuple of integers' + ) + + # ===================================================================================================================== + # Define the description + # ===================================================================================================================== + @property + def description(self): + return self._description + + @description.setter + def description(self, val): + if isinstance(val, str): + self._description = val + else: + raise ValueError('Invalid description. Must be a string.') + + +class TypeCheckedList(list): + def __init__(self, checkItem, itemList=None): + super(TypeCheckedList, self).__init__() + self.checkItem = checkItem + + if itemList is not None: + if isinstance(itemList, list) or isinstance(itemList, tuple): + for itm in itemList: + self.append(itm) + else: + raise ValueError('Input to itemList is not iterable') + + def __setitem__(self, key, val): + if isinstance(val, self.checkItem): + super(TypeCheckedList, self).__setitem__(key, val) + elif isinstance(val, (tuple, list)): + cks = [isinstance(vl, self.checkItem) for vl in val] + if sum(cks) == len(cks): + super(TypeCheckedList, self).__setitem__(key, val) + else: + raise ValueError('Input must be an instance of the defined type') + else: + raise ValueError('Input must be an instance of the defined type') + + def append(self, val): + if isinstance(val, self.checkItem): + super(TypeCheckedList, self).append(val) + else: + raise ValueError('Input must be an instance of the defined type') + + +class BBList(TypeCheckedList): + def __init__(self): + super(BBList, self).__init__(BlackBoxFunctionModel_Variable, []) + self._lookupDict = {} + self._counter = 0 + + def __getitem__(self, val): + if isinstance(val, int): + return super(BBList, self).__getitem__(val) + elif isinstance(val, str): + return super(BBList, self).__getitem__(self._lookupDict[val]) + else: + raise ValueError('Input must be an integer or a valid variable name') + + def append(*args, **kwargs): + args = list(args) + self = args.pop(0) + + if len(args) + len(kwargs.values()) == 1: + if len(args) == 1: + inputData = args[0] + if len(kwargs.values()) == 1: + inputData = list(kwargs.values())[0] + + if isinstance(inputData, self.checkItem): + if inputData.name in self._lookupDict.keys(): + raise ValueError( + "Key '%s' already exists in the input list" % (inputData.name) + ) + self._lookupDict[inputData.name] = self._counter + self._counter += 1 + super(BBList, self).append(inputData) + else: + if isinstance(inputData, str): + raise ValueError( + "Key '%s' not passed in to the black box variable constructor" + % ('units') + ) + else: + raise ValueError('Invalid (single) input type') + + elif len(args) + len(kwargs.values()) <= 4: + argKeys = ['name', 'units', 'description', 'size'] + ipd = dict(zip(argKeys[0 : len(args)], args)) + for ky, vl in kwargs.items(): + if ky in ipd: + raise ValueError( + "Key '%s' declared after non-keyword arguments and is out of order" + % (ky) + ) + else: + ipd[ky] = vl + + for ak in argKeys: + if ak not in ipd.keys(): + if ak == 'description': + ipd['description'] = '' + elif ak == 'size': + ipd['size'] = 0 + else: + raise ValueError( + "Key '%s' not passed in to the black box variable constructor" + % (ak) + ) + + if ipd['name'] in self._lookupDict.keys(): + raise ValueError( + "Key '%s' already exists in the input list" % (ipd['name']) + ) + self._lookupDict[ipd['name']] = self._counter + self._counter += 1 + super(BBList, self).append(BlackBoxFunctionModel_Variable(**ipd)) + + else: + raise ValueError('Too many inputs to a black box variable') + + +errorString = 'This function is calling to the base class and has not been defined.' + + +def toList(itm, keylist): + if isinstance(itm, (tuple, list)): + return list(itm) + elif isinstance(itm, dict): + retList = [] + for ky in keylist: + if ky in itm.keys(): + retList.append(itm[ky]) + return retList + else: + return [itm] + + +class BlackBoxFunctionModel(ExternalGreyBoxModel): + # --------------------------------------------------------------------------------------------------------------------- + # --------------------------------------------------------------------------------------------------------------------- + def __init__(self): + # pyomo.common.numeric_types.RegisterNumericType(pint.Quantity) + super(BlackBoxFunctionModel, self).__init__() + + # List of the inputs and outputs + self.inputs = BBList() + self.outputs = BBList() + + self.inputVariables_optimization = None + self.outputVariables_optimization = None + + # A simple description of the model + self.description = None + + # Defines the order of derivative available in the black box + self.availableDerivative = 0 + + self._cache = None + self._NunwrappedOutputs = None + self._NunwrappedInputs = None + + def setOptimizationVariables( + self, inputVariables_optimization, outputVariables_optimization + ): + self.inputVariables_optimization = inputVariables_optimization + self.outputVariables_optimization = outputVariables_optimization + + # --------------------------------------------------------------------------------------------------------------------- + # pyomo things + # --------------------------------------------------------------------------------------------------------------------- + def input_names(self): + inputs_unwrapped = [] + for ivar in self.inputVariables_optimization: + if isinstance(ivar, pyomo.core.base.var.ScalarVar): + inputs_unwrapped.append(ivar) + elif isinstance(ivar, pyomo.core.base.var.IndexedVar): + validIndices = list(ivar.index_set().data()) + for vi in validIndices: + inputs_unwrapped.append(ivar[vi]) + else: + raise ValueError("Invalid type for input variable") + + return [ip.__str__() for ip in inputs_unwrapped] + + def output_names(self): + outputs_unwrapped = [] + for ovar in self.outputVariables_optimization: + if isinstance(ovar, pyomo.core.base.var.ScalarVar): + outputs_unwrapped.append(ovar) + elif isinstance(ovar, pyomo.core.base.var.IndexedVar): + validIndices = list(ovar.index_set().data()) + for vi in validIndices: + outputs_unwrapped.append(ovar[vi]) + else: + raise ValueError("Invalid type for output variable") + + return [op.__str__() for op in outputs_unwrapped] + + def set_input_values(self, input_values): + self._input_values = input_values + self._cache = None + + def evaluate_outputs(self): + self.fillCache() + opts = self._cache['pyomo_outputs'] + return opts + + def evaluate_jacobian_outputs(self): + self.fillCache() + jac = self._cache['pyomo_jacobian'] + return jac + + def post_init_setup(self, defaultVal=1.0): + self._input_values = np.ones(self._NunwrappedInputs) * defaultVal + + def fillCache(self): + if self._cache is None: + self._cache = {} + + raw_inputs = self._input_values + bb_inputs = [] + + ptr = 0 + + for i in range(0, len(self.inputVariables_optimization)): + optimizationInput = self.inputVariables_optimization[i] + if not isinstance( + optimizationInput, + (pyomo.core.base.var.IndexedVar, pyomo.core.base.var.ScalarVar), + ): + raise ValueError("Invalid input variable type") + + ipt = self.inputs[i] + + shape = [len(idx) for idx in optimizationInput.index_set().subsets()] + localShape = ipt.size + + optimizationUnits = self.inputVariables_optimization[i].get_units() + localUnits = ipt.units + + if isinstance(optimizationInput, pyomo.core.base.var.IndexedVar): + value = np.zeros(shape) + for vix in list(optimizationInput.index_set().data()): + raw_val = float(raw_inputs[ptr]) * optimizationUnits + raw_val_correctedUnits = pyomo_units.convert( + raw_val, localUnits + ) + value[vix] = pyo.value(raw_val_correctedUnits) + ptr += 1 + self.sizeCheck(localShape, value * localUnits) + bb_inputs.append(value * localUnits) + + else: # isinstance(optimizationInput, pyomo.core.base.var.ScalarVar):, just checked this + value = raw_inputs[ptr] * optimizationUnits + ptr += 1 + self.sizeCheck(localShape, value) + value_correctedUnits = pyomo_units.convert(value, localUnits) + bb_inputs.append(value_correctedUnits) + + bbo = self.BlackBox(*bb_inputs) + + self._cache['raw'] = bbo + self._cache['raw_value'] = bbo[0] + self._cache['raw_jacobian'] = bbo[1] + + outputVector = [] + if not isinstance(bbo[0], (list, tuple)): + valueList = [bbo[0]] + jacobianList = [bbo[1]] + else: + valueList = bbo[0] + jacobianList = bbo[1] + for i in range(0, len(valueList)): + optimizationOutput = self.outputVariables_optimization[i] + if not isinstance( + optimizationOutput, + (pyomo.core.base.var.IndexedVar, pyomo.core.base.var.ScalarVar), + ): + raise ValueError("Invalid output variable type") + opt = self.outputs[i] + + modelOutputUnits = opt.units + outputOptimizationUnits = optimizationOutput.get_units() + vl = valueList[i] + if isinstance(vl, NumericNDArray): + validIndexList = optimizationOutput.index_set().data() + for j in range(0, len(validIndexList)): + vi = validIndexList[j] + corrected_value = pyo.value( + pyomo_units.convert(vl[vi], outputOptimizationUnits) + ) # now unitless in correct units + outputVector.append(corrected_value) + + elif isinstance( + vl, + ( + pyomo.core.expr.numeric_expr.NPV_ProductExpression, + pyomo.core.base.units_container._PyomoUnit, + ), + ): + corrected_value = pyo.value( + pyomo_units.convert(vl, outputOptimizationUnits) + ) # now unitless in correct units + outputVector.append(corrected_value) + + else: + raise ValueError("Invalid output variable type") + + self._cache['pyomo_outputs'] = outputVector + + outputJacobian = ( + np.ones([self._NunwrappedOutputs, self._NunwrappedInputs]) * -1 + ) + ptr_row = 0 + ptr_col = 0 + + for i in range(0, len(jacobianList)): + oopt = self.outputVariables_optimization[i] + # Checked about 20 lines above + # if not isinstance(oopt, (pyomo.core.base.var.ScalarVar,pyomo.core.base.var.IndexedVar)): + # raise ValueError("Invalid type for output variable") + lopt = self.outputs[i] + oounits = oopt.get_units() + lounits = lopt.units + # oshape = [len(idx) for idx in oopt.index_set().subsets()] + ptr_col = 0 + for j in range(0, len(self.inputs)): + oipt = self.inputVariables_optimization[j] + # This is checked about 80 lines up + # if not isinstance(oipt, (pyomo.core.base.var.ScalarVar,pyomo.core.base.var.IndexedVar)): + # raise ValueError("Invalid type for output variable") + lipt = self.inputs[j] + oiunits = oipt.get_units() + liunits = lipt.units + # ishape = [len(idx) for idx in oipt.index_set().subsets()] + + jacobianValue_raw = jacobianList[i][j] + + if isinstance( + jacobianValue_raw, + ( + pyomo.core.expr.numeric_expr.NPV_ProductExpression, + pyomo.core.base.units_container._PyomoUnit, + ), + ): + corrected_value = pyo.value( + pyomo_units.convert(jacobianValue_raw, oounits / oiunits) + ) # now unitless in correct units + outputJacobian[ptr_row, ptr_col] = corrected_value + ptr_col += 1 + ptr_row_step = 1 + + elif isinstance(jacobianValue_raw, NumericNDArray): + jshape = jacobianValue_raw.shape + + if isinstance(oopt, pyomo.core.base.var.ScalarVar): + oshape = 0 + else: # isinstance(oopt, pyomo.core.base.var.IndexedVar), checked above + oshape = [len(idx) for idx in oopt.index_set().subsets()] + + if isinstance(oipt, pyomo.core.base.var.ScalarVar): + ishape = 0 + else: # isinstance(oipt, pyomo.core.base.var.IndexedVar), checked above + ishape = [len(idx) for idx in oipt.index_set().subsets()] + + if oshape == 0: + validIndices = list(oipt.index_set().data()) + for vi in validIndices: + corrected_value = pyo.value( + pyomo_units.convert( + jacobianValue_raw[vi], oounits / oiunits + ) + ) # now unitless in correct units + outputJacobian[ptr_row, ptr_col] = corrected_value + ptr_col += 1 + ptr_row_step = 1 + + elif ishape == 0: + ptr_row_cache = ptr_row + validIndices = list(oopt.index_set().data()) + for vi in validIndices: + corrected_value = pyo.value( + pyomo_units.convert( + jacobianValue_raw[vi], oounits / oiunits + ) + ) # now unitless in correct units + outputJacobian[ptr_row, ptr_col] = corrected_value + ptr_row += 1 + ptr_row = ptr_row_cache + ptr_row_step = len(validIndices) + + # elif ishape == 0 and oshape == 0: # Handled by the scalar case above + + else: + # both are dimensioned vectors + # oshape, ishape, jshape + ptr_row_cache = ptr_row + ptr_col_cache = ptr_col + validIndices_o = list(oopt.index_set().data()) + validIndices_i = list(oipt.index_set().data()) + + for vio in validIndices_o: + if isinstance(vio, (float, int)): + vio = (vio,) + for vii in validIndices_i: + if isinstance(vii, (float, int)): + vii = (vii,) + corrected_value = pyo.value( + pyomo_units.convert( + jacobianValue_raw[vio + vii], + oounits / oiunits, + ) + ) # now unitless in correct units + outputJacobian[ptr_row, ptr_col] = corrected_value + ptr_col += 1 + ptr_col = ptr_col_cache + ptr_row += 1 + ptr_row = ptr_row_cache + ptr_row_step = len(validIndices_o) + + else: + raise ValueError("Invalid jacobian type") + + ptr_row += ptr_row_step + + # print(outputJacobian.dtype) + # print(sps.coo_matrix(outputJacobian)) + # outputJacobian = np.zeros([3,1]) + self._cache['pyomo_jacobian'] = sps.coo_matrix( + outputJacobian, shape=outputJacobian.shape + ) + + # --------------------------------------------------------------------------------------------------------------------- + # --------------------------------------------------------------------------------------------------------------------- + # These models must be defined in each individual model, just placeholders here + def BlackBox(*args, **kwargs): + raise AttributeError(errorString) + + # --------------------------------------------------------------------------------------------------------------------- + # --------------------------------------------------------------------------------------------------------------------- + # These functions wrap black box to provide more functionality + + def BlackBox_Standardized(self, *args, **kwargs): + runCases, extras = self.parseInputs(*args, **kwargs) + if len(runCases) != 1: + # This is actually a MultiCase, run the correct function instead + return self.MultiCase(*args, **kwargs) + + # -------- + # Run + # -------- + opt_raw = self.BlackBox(*args, **kwargs) + + # -------- + # Setup + # -------- + outputNames = [o.name for o in self.outputs] + inputNames = [i.name for i in self.inputs] + + returnTuple = namedtuple('returnTuple', ['values', 'first', 'second']) + optTuple = namedtuple('optTuple', outputNames) + iptTuple = namedtuple('iptTuple', inputNames) + + structuredOutput = {} + structuredOutput['values'] = np.zeros(len(self.outputs)).tolist() + structuredOutput['first'] = np.zeros( + [len(self.outputs), len(self.inputs)] + ).tolist() + structuredOutput['second'] = np.zeros( + [len(self.outputs), len(self.inputs), len(self.inputs)] + ).tolist() + + opt = toList(opt_raw, ['values', 'first', 'second']) + # print(opt) + + # -------- + # Values + # -------- + # values will always come back as an appropriate list + opt_values = toList(opt[0], []) + for i in range(0, len(structuredOutput['values'])): + if self.outputs[i].size == 0: + structuredOutput['values'][i] = 0.0 * self.outputs[i].units + else: + currentSize = toList(self.outputs[i].size, []) + structuredOutput['values'][i] = ( + np.zeros(currentSize) * self.outputs[i].units + ) + + for i, vl in enumerate(opt_values): + if self.outputs[i].size == 0: + structuredOutput['values'][i] = self.convert( + opt_values[i], self.outputs[i].units + ) + elif isinstance(self.outputs[i].size, int): + for ix in range(0, self.outputs[i].size): + structuredOutput['values'][i][ix] = opt_values[i][ix] + else: + listOfIndices = list( + itertools.product(*[range(0, n) for n in self.outputs[i].size]) + ) + for ix in listOfIndices: + structuredOutput['values'][i][ix] = opt_values[i][ + ix + ] # unit conversion handled automatically by pint + + # structuredOutput['values'] = opt_values + + # -------- + # First + # -------- + # first orders will be [num] or [num, num, lst, num...] or [lst, lst, lst, lst...] + # should be [lst, lst, lst, lst...] + if len(opt) > 1: + opt_first = toList(opt[1], outputNames) + if len(self.outputs) == 1 and len(opt_first) == len(self.inputs): + opt_first = [opt_first] + for i, vl in enumerate(opt_first): + opt_first[i] = toList(vl, inputNames) + + for i in range(0, len(structuredOutput['first'])): + currentSize_output = toList(self.outputs[i].size, []) + for j in range(0, len(structuredOutput['first'][i])): + if self.outputs[i].size == 0 and self.inputs[j].size == 0: + structuredOutput['first'][i][j] = ( + 0.0 * self.outputs[i].units / self.inputs[j].units + ) + elif self.outputs[i].size == 0: + currentSize_input = toList(self.inputs[j].size, []) + structuredOutput['first'][i][j] = ( + np.zeros(currentSize_input) + * self.outputs[i].units + / self.inputs[j].units + ) + elif self.inputs[j].size == 0: + structuredOutput['first'][i][j] = ( + np.zeros(currentSize_output) + * self.outputs[i].units + / self.inputs[j].units + ) + else: + currentSize_input = toList(self.inputs[j].size, []) + structuredOutput['first'][i][j] = ( + np.zeros(currentSize_output + currentSize_input) + * self.outputs[i].units + / self.inputs[j].units + ) + + for i, vl in enumerate(opt_first): + for j, vlj in enumerate(opt_first[i]): + # from pyomo.contrib.edi import ediprint + # print() + # ediprint(opt_first) + # ediprint(opt_first[i]) + # ediprint(opt_first[i][j]) + # print(self.outputs) + # print(self.inputs) + # print(i) + # print(j) + if self.outputs[i].size == 0 and self.inputs[j].size == 0: + structuredOutput['first'][i][j] = self.convert( + opt_first[i][j], + self.outputs[i].units / self.inputs[j].units, + ) + elif self.outputs[i].size == 0: + # print(self.inputs[j].size) + if isinstance(self.inputs[j].size, int): + sizelist = [self.inputs[j].size] + else: + assert isinstance(self.inputs[j].size, list) + sizelist = self.inputs[j].size + listOfIndices = list( + itertools.product(*[range(0, n) for n in sizelist]) + ) + # print(listOfIndices) + for ix in listOfIndices: + if len(ix) == 1: + ix = ix[0] + structuredOutput['first'][i][j][ix] = opt_first[i][j][ + ix + ] # unit conversion handled automatically by pint + elif self.inputs[j].size == 0: + if isinstance(self.outputs[i].size, int): + sizelist = [self.outputs[i].size] + else: + assert isinstance(self.outputs[i].size, list) + sizelist = self.outputs[i].size + listOfIndices = list( + itertools.product(*[range(0, n) for n in sizelist]) + ) + for ix in listOfIndices: + if len(ix) == 1: + ix = ix[0] + structuredOutput['first'][i][j][ix] = opt_first[i][j][ + ix + ] # unit conversion handled automatically by pint + else: + if isinstance(self.inputs[j].size, int): + sizelist_inputs = [self.inputs[j].size] + else: + assert isinstance(self.inputs[j].size, list) + sizelist_inputs = self.inputs[j].size + + if isinstance(self.outputs[i].size, int): + sizelist_outputs = [self.outputs[i].size] + else: + assert isinstance(self.inputs[j].size, list) + sizelist_outputs = self.outputs[i].size + + listOfIndices = list( + itertools.product( + *[ + range(0, n) + for n in sizelist_outputs + sizelist_inputs + ] + ) + ) + for ix in listOfIndices: + if len(ix) == 1: + ix = ix[0] + structuredOutput['first'][i][j][ix] = opt_first[i][j][ + ix + ] # unit conversion handled automatically by pint + else: + structuredOutput['first'] = None + # for i,vl in enumerate(opt_first): + # opt_first[i] = toList(vl,inputNames) + # structuredOutput['first'] = opt_first + + # -------- + # Second + # -------- + # if len(opt)>2: + # opt_second = toList(opt[2],outputNames) + # # Not implemented in this version + # else: + structuredOutput['second'] = None + + # for i,vl in enumerate(opt_second): + # opt_second[i] = toList(vl,inputNames) + # for j, vlj in enumerate(opt_second[i]): + # opt_second[i][j] = toList(vlj,inputNames) + # structuredOutput['second'] = opt_second + + # -------- + # Pack + # -------- + # for i in range(0, len(structuredOutput['second'])): + # for j in range(0, len(structuredOutput['second'][i])): + # structuredOutput['second'][i][j] = iptTuple( + # *structuredOutput['second'][i][j] + # ) + + for i in range(0, len(structuredOutput['first'])): + structuredOutput['first'][i] = iptTuple(*structuredOutput['first'][i]) + # structuredOutput['second'][i] = iptTuple(*structuredOutput['second'][i]) + + structuredOutput['values'] = optTuple(*structuredOutput['values']) + structuredOutput['first'] = optTuple(*structuredOutput['first']) + # structuredOutput['second'] = optTuple(*structuredOutput['second']) + + return returnTuple( + structuredOutput['values'], + structuredOutput['first'], + structuredOutput['second'], + ) + + def MultiCase(self, *args, **kwargs): + runCases, extras = self.parseInputs(*args, **kwargs) + + outputList = [] + # values = [] + # first = [] + # second = [] + for rc in runCases: + try: + ipt = rc | extras + except: + # Allow for python 3.8 compatibility + ipt = copy.deepcopy(rc) + ipt.update(extras) + + opt = self.BlackBox_Standardized(**ipt) + outputList.append(opt) + # values.append(opt[0]) + # first.append(opt[1]) + # second.append(opt[2]) + + # returnTuple = namedtuple('returnTuple', ['values','first','second']) + # return returnTuple(values, first, second) + return outputList + + # --------------------------------------------------------------------------------------------------------------------- + # --------------------------------------------------------------------------------------------------------------------- + # unit handling functions + + def convert(self, val, unts): + if isinstance(val, pint.Quantity): + val.to(pyomo_units._get_pint_units(unts)) + return val + + try: + val = val * pyomo_units.dimensionless + val = pyo.value(val) * pyomo_units.get_units(val) + except: + pass ## will handle later + + if isinstance( + val, + ( + pyomo.core.base.units_container._PyomoUnit, + pyomo.core.expr.numeric_expr.NPV_ProductExpression, + ), + ): + rval = pyomo_units.convert(val, unts) + return pyo.value(rval) * pyomo_units.get_units(rval) + elif isinstance(val, NumericNDArray): + shp = val.shape + ix = np.ndindex(*shp) + opt = np.zeros(shp) + for i in range(0, np.prod(shp)): + ixt = next(ix) + opt[ixt] = pyo.value(pyomo_units.convert(val[ixt], unts)) + return opt * unts + # elif isinstance(val, (list, tuple)): + else: + raise ValueError('Invalid type passed to unit conversion function') + + def pyomo_value(self, val): + try: + return pyo.value(val) + except: + if isinstance(val, NumericNDArray): + shp = val.shape + ix = np.ndindex(*shp) + opt = np.zeros(shp) + for i in range(0, np.prod(shp)): + ixt = next(ix) + opt[ixt] = pyo.value(val[ixt]) + return opt + else: + raise ValueError('Invalid type passed to pyomo_value function') + + # --------------------------------------------------------------------------------------------------------------------- + # --------------------------------------------------------------------------------------------------------------------- + def parseInputs(self, *args, **kwargs): + args = list(args) # convert tuple to list + inputNames = [self.inputs[i].name for i in range(0, len(self.inputs))] + + # ------------------------------ + # ------------------------------ + if len(args) + len(kwargs.values()) == 1: + if len(args) == 1: + inputData = args[0] + if len(kwargs.values()) == 1: + inputData = list(kwargs.values())[0] + + if len(inputNames) == 1: + try: + rs = self.sanitizeInputs(inputData) + return ( + [dict(zip(inputNames, [rs]))], + # -self.availableDerivative - 1, + {}, + ) # one input being passed in + except: + pass # otherwise, proceed + + if isinstance(inputData, (list, tuple)): + dataRuns = [] + for idc in inputData: + if isinstance(idc, dict): + sips = self.sanitizeInputs(**idc) + if len(inputNames) == 1: + sips = [sips] + runDictS = dict(zip(inputNames, sips)) + dataRuns.append( + runDictS + ) # the BlackBox([{'x1':x1, 'x2':x2},{'x1':x1, 'x2':x2},...]) case + elif isinstance(idc, (list, tuple)): + if len(idc) == len(inputNames): + sips = self.sanitizeInputs(*idc) + if len(inputNames) == 1: + sips = [sips] + runDictS = dict(zip(inputNames, sips)) + dataRuns.append( + runDictS + ) # the BlackBox([ [x1, x2], [x1, x2],...]) case + else: + raise ValueError( + 'Entry in input data list has improper length' + ) + else: + raise ValueError( + "Invalid data type in the input list. Note that BlackBox([x1,x2,y]) must be passed in as BlackBox([[x1,x2,y]]) or " + + "BlackBox(*[x1,x2,y]) or BlackBox({'x1':x1,'x2':x2,'y':y}) or simply BlackBox(x1, x2, y) to avoid processing singularities. " + + "Best practice is BlackBox({'x1':x1,'x2':x2,'y':y}). For multi-case input, be sure to wrap even single inputs in lists," + + ' ex: [[x],[x],[x]] and not [x,x,x].' + ) + return dataRuns, {} + + elif isinstance(inputData, dict): + if set(list(inputData.keys())) == set(inputNames): + try: + inputLengths = [len(inputData[kw]) for kw in inputData.keys()] + except: + sips = self.sanitizeInputs(**inputData) + if len(inputNames) == 1: + sips = [sips] + return ( + [dict(zip(inputNames, sips))], + # -self.availableDerivative - 1, + {}, + ) # the BlackBox(*{'x1':x1, 'x2':x2}) case, somewhat likely + + if not all( + [ + inputLengths[i] == inputLengths[0] + for i in range(0, len(inputLengths)) + ] + ): + sips = self.sanitizeInputs(**inputData) + return ( + [dict(zip(inputNames, sips))], + # -self.availableDerivative - 1, + {}, + ) # the BlackBox(*{'x1':x1, 'x2':x2}) case where vectors for x1... are passed in (likely to fail on previous line) + else: + try: + sips = self.sanitizeInputs(**inputData) + return ( + [dict(zip(inputNames, sips))], + # -self.availableDerivative - 1, + {}, + ) # the BlackBox(*{'x1':x1, 'x2':x2}) case where vectors all inputs have same length intentionally (likely to fail on previous line) + except: + dataRuns = [] + for i in range(0, inputLengths[0]): + runDict = {} + for ky, vl in inputData.items(): + runDict[ky] = vl[i] + sips = self.sanitizeInputs(**runDict) + if len(inputNames) == 1: + sips = [sips] + runDictS = dict(zip(inputNames, sips)) + dataRuns.append(runDictS) + return ( + dataRuns, + # self.availableDerivative, + {}, + ) # the BlackBox({'x1':x1_vec, 'x2':x2_vec}) case, most likely + else: + raise ValueError('Keywords did not match the expected list') + else: + raise ValueError('Got unexpected data type %s' % (str(type(inputData)))) + # ------------------------------ + # ------------------------------ + else: + if any( + [ + list(kwargs.keys())[i] in inputNames + for i in range(0, len(list(kwargs.keys()))) + ] + ): + # some of the inputs are defined in the kwargs + if len(args) >= len(inputNames): + raise ValueError( + 'A keyword input is defining an input, but there are too many unkeyed arguments for this to occur. Check the inputs.' + ) + else: + if len(args) != 0: + availableKeywords = inputNames[-len(args) :] + else: + availableKeywords = inputNames + + valList = args + [None] * (len(inputNames) - len(args)) + for ky in availableKeywords: + ix = inputNames.index(ky) + valList[ix] = kwargs[ky] + + # Not possible to reach due to other checks + # if any([valList[i]==None for i in range(0,len(valList))]): + # raise ValueError('Keywords did not properly fill in the remaining arguments. Check the inputs.') + + sips = self.sanitizeInputs(*valList) + + if len(inputNames) == 1: + sips = [sips] + + remainingKwargs = copy.deepcopy(kwargs) + for nm in inputNames: + try: + del remainingKwargs[nm] + except: + # was in args + pass + + return ( + [dict(zip(inputNames, sips))], + # -self.availableDerivative - 1, + remainingKwargs, + ) # Mix of args and kwargs define inputs + else: + # all of the inputs are in the args + try: + sips = self.sanitizeInputs(*args[0 : len(inputNames)]) + + if len(inputNames) == 1: + sips = [sips] + + remainingKwargs = copy.deepcopy(kwargs) + remainingKwargs['remainingArgs'] = args[len(inputNames) :] + return ( + [dict(zip(inputNames, sips))], + # -self.availableDerivative - 1, + remainingKwargs, + ) # all inputs are in args + # except: + except Exception as e: + # not possible due to other checks + # if str(e) == 'Too many inputs': + # raise ValueError(e) + if str(e) == 'Not enough inputs': + raise ValueError(e) + else: # otherwise, proceed + runCases, extra_singleInput = self.parseInputs(args[0]) + remainingKwargs = copy.deepcopy(kwargs) + remainingKwargs['remainingArgs'] = args[len(inputNames) :] + return runCases, remainingKwargs + + # --------------------------------------------------------------------------------------------------------------------- + # --------------------------------------------------------------------------------------------------------------------- + def sizeCheck(self, size, ipval_correctUnits): + if size is not None: + szVal = ipval_correctUnits + if isinstance( + szVal, + ( + pyomo.core.expr.numeric_expr.NPV_ProductExpression, + pyomo.core.base.units_container._PyomoUnit, + ), + ): + if size != 0 and size != 1: + raise ValueError( + 'Size did not match the expected size %s (ie: Scalar)' + % (str(size)) + ) + elif isinstance(szVal, NumericNDArray): + shp = szVal.shape + if isinstance(size, (int, float)): + size = [size] + # else: + if len(shp) != len(size): + raise ValueError( + 'Shapes/Sizes of %s does not match the expected %s' + % (str(shp), str(size)) + ) + for j in range(0, len(shp)): + if size[j] != -1: # was declared of flexible length + if size[j] != shp[j]: + raise ValueError( + 'Shapes/Sizes of %s does not match the expected %s' + % (str(shp), str(size)) + ) + else: + raise ValueError('Invalid type detected when checking size') + + # --------------------------------------------------------------------------------------------------------------------- + # --------------------------------------------------------------------------------------------------------------------- + def sanitizeInputs(self, *args, **kwargs): + nameList = [self.inputs[i].name for i in range(0, len(self.inputs))] + + if len(args) + len(kwargs.values()) > len(nameList): + raise ValueError('Too many inputs') + if len(args) + len(kwargs.values()) < len(nameList): + raise ValueError('Not enough inputs') + inputDict = {} + for i in range(0, len(args)): + rg = args[i] + inputDict[nameList[i]] = rg + + for ky, vl in kwargs.items(): + if ky in nameList: + inputDict[ky] = vl + else: + raise ValueError( + 'Unexpected input keyword argument %s in the inputs' % (ky) + ) + + opts = [] + + for i in range(0, len(nameList)): + name = nameList[i] + nameCheck = self.inputs[i].name + unts = self.inputs[i].units + size = self.inputs[i].size + + # should be impossible + # if name != nameCheck: + # raise RuntimeError('Something went wrong and values are not consistent. Check your defined inputs.') + + ipval = inputDict[name] + + if isinstance(ipval, NumericNDArray): + for ii in range(0, len(ipval)): + try: + ipval[ii] = self.convert(ipval[ii], unts) # ipval.to(unts) + except: + raise ValueError( + 'Could not convert %s of %s to %s' + % (name, str(ipval), str(unts)) + ) + ipval_correctUnits = ipval + else: + try: + ipval_correctUnits = self.convert(ipval, unts) # ipval.to(unts) + except: + raise ValueError( + 'Could not convert %s of %s to %s' + % (name, str(ipval), str(unts)) + ) + + # superseded by the custom convert function + # if not isinstance(ipval_correctUnits, (pyomo.core.expr.numeric_expr.NPV_ProductExpression, + # NumericNDArray, + # pyomo.core.base.units_container._PyomoUnit)): + # ipval_correctUnits = ipval_correctUnits * pyomo_units.dimensionless + + self.sizeCheck(size, ipval_correctUnits) + + opts.append(ipval_correctUnits) + if len(opts) == 1: + opts = opts[0] + + return opts + + # --------------------------------------------------------------------------------------------------------------------- + # --------------------------------------------------------------------------------------------------------------------- + def checkOutputs(self, *args, **kwargs): + raise NotImplementedError('Contact developers to use this function') + # nameList = [self.outputs[i].name for i in range(0,len(self.outputs))] + # if len(args) + len(kwargs.values()) > len(nameList): + # raise ValueError('Too many outputs') + # if len(args) + len(kwargs.values()) < len(nameList): + # raise ValueError('Not enough outputs') + # inputDict = {} + # for i in range(0,len(args)): + # rg = args[i] + # inputDict[nameList[i]] = rg + + # for ky, vl in kwargs.items(): + # if ky in nameList: + # inputDict[ky] = vl + # else: + # raise ValueError('Unexpected output keyword argument %s in the outputs'%(ky)) + + # opts = [] + + # for i in range(0,len(nameList)): + # name = nameList[i] + # nameCheck = self.outputs[i].name + # unts = self.outputs[i].units + # size = self.outputs[i].size + + # if name != nameCheck: + # raise RuntimeError('Something went wrong and values are not consistent. Check your defined inputs.') + + # ipval = inputDict[name] + + # if unts is not None: + # try: + # ipval_correctUnits = self.convert(ipval, unts) + # except: + # raise ValueError('Could not convert %s of %s to %s'%(name, str(ipval),str(unts))) + # else: + # ipval_correctUnits = ipval + + # if not isinstance(ipval_correctUnits, (pyomo.core.expr.numeric_expr.NPV_ProductExpression, + # NumericNDArray, + # pyomo.core.base.units_container._PyomoUnit)): + # ipval_correctUnits = ipval_correctUnits * pyomo_units.dimensionless + + # self.sizeCheck(size, ipval_correctUnits) + + # opts.append(ipval_correctUnits) + # if len(opts) == 1: + # opts = opts[0] + + # return opts + + # --------------------------------------------------------------------------------------------------------------------- + # --------------------------------------------------------------------------------------------------------------------- + def getSummary(self, whitespace=6): + pstr = '\n' + pstr += 'Model Description\n' + pstr += '=================\n' + descr_str = self.description.__repr__() + pstr += descr_str[1:-1] + '\n\n' + + longestName = 0 + longestUnits = 0 + longestSize = 0 + for ipt in self.inputs: + nml = len(ipt.name) + if nml > longestName: + longestName = nml + unts = ipt.units.__str__() # _repr_html_() + # unts = unts.replace('','^') + # unts = unts.replace('','') + # unts = unts.replace('\[', '[') + # unts = unts.replace('\]', ']') + unl = len(unts) + if unl > longestUnits: + longestUnits = unl + if type(ipt.size) == list: + lsz = len(ipt.size.__repr__()) + else: + lsz = len(str(ipt.size)) + if lsz > longestSize: + longestSize = lsz + namespace = max([4, longestName]) + whitespace + unitspace = max([5, longestUnits]) + whitespace + sizespace = max([4, longestSize]) + whitespace + fulllength = namespace + unitspace + sizespace + 11 + pstr += 'Inputs\n' + pstr += '=' * fulllength + pstr += '\n' + pstr += 'Name'.ljust(namespace) + pstr += 'Units'.ljust(unitspace) + pstr += 'Size'.ljust(sizespace) + pstr += 'Description' + pstr += '\n' + pstr += '-' * (namespace - whitespace) + pstr += ' ' * whitespace + pstr += '-' * (unitspace - whitespace) + pstr += ' ' * whitespace + pstr += '-' * (sizespace - whitespace) + pstr += ' ' * whitespace + pstr += '-----------' + pstr += '\n' + for ipt in self.inputs: + pstr += ipt.name.ljust(namespace) + unts = ipt.units.__str__() # _repr_html_() + # unts = unts.replace('','^') + # unts = unts.replace('','') + # unts = unts.replace('\[', '[') + # unts = unts.replace('\]', ']') + pstr += unts.ljust(unitspace) + lnstr = '%s' % (ipt.size.__repr__()) + pstr += lnstr.ljust(sizespace) + pstr += ipt.description + pstr += '\n' + pstr += '\n' + + longestName = 0 + longestUnits = 0 + longestSize = 0 + for opt in self.outputs: + nml = len(opt.name) + if nml > longestName: + longestName = nml + unts = opt.units.__str__() # _repr_html_() + # unts = unts.replace('','^') + # unts = unts.replace('','') + # unts = unts.replace('\[', '[') + # unts = unts.replace('\]', ']') + unl = len(unts) + if unl > longestUnits: + longestUnits = unl + if type(opt.size) == list: + lsz = len(opt.size.__repr__()) + else: + lsz = len(str(opt.size)) + if lsz > longestSize: + longestSize = lsz + namespace = max([4, longestName]) + whitespace + unitspace = max([5, longestUnits]) + whitespace + sizespace = max([4, longestSize]) + whitespace + fulllength = namespace + unitspace + sizespace + 11 + pstr += 'Outputs\n' + pstr += '=' * fulllength + pstr += '\n' + pstr += 'Name'.ljust(namespace) + pstr += 'Units'.ljust(unitspace) + pstr += 'Size'.ljust(sizespace) + pstr += 'Description' + pstr += '\n' + pstr += '-' * (namespace - whitespace) + pstr += ' ' * whitespace + pstr += '-' * (unitspace - whitespace) + pstr += ' ' * whitespace + pstr += '-' * (sizespace - whitespace) + pstr += ' ' * whitespace + pstr += '-----------' + pstr += '\n' + for opt in self.outputs: + pstr += opt.name.ljust(namespace) + unts = opt.units.__str__() # _repr_html_() + # unts = unts.replace('','^') + # unts = unts.replace('','') + # unts = unts.replace('\[', '[') + # unts = unts.replace('\]', ']') + pstr += unts.ljust(unitspace) + lnstr = '%s' % (opt.size.__repr__()) + pstr += lnstr.ljust(sizespace) + pstr += opt.description + pstr += '\n' + pstr += '\n' + + return pstr + + # --------------------------------------------------------------------------------------------------------------------- + # --------------------------------------------------------------------------------------------------------------------- + + @property + def summary(self): + return self.getSummary() + + # --------------------------------------------------------------------------------------------------------------------- + # --------------------------------------------------------------------------------------------------------------------- + def __repr__(self): + pstr = 'AnalysisModel( [' + for i in range(0, len(self.outputs)): + pstr += self.outputs[i].name + pstr += ',' + pstr = pstr[0:-1] + pstr += ']' + pstr += ' == ' + pstr += 'f(' + for ipt in self.inputs: + pstr += ipt.name + pstr += ', ' + pstr = pstr[0:-2] + pstr += ')' + pstr += ' , ' + pstr = pstr[0:-2] + pstr += '])' + return pstr diff --git a/pyomo/contrib/edi/examples/__init__.py b/pyomo/contrib/edi/examples/__init__.py new file mode 100644 index 00000000000..38b30839be3 --- /dev/null +++ b/pyomo/contrib/edi/examples/__init__.py @@ -0,0 +1,16 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# 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. +# +# Development of this module was conducted as part of the Institute for +# the Design of Advanced Energy Systems (IDAES) with support through the +# Simulation-Based Engineering, Crosscutting Research Program within the +# U.S. Department of Energy’s Office of Fossil Energy and Carbon Management. +# +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ diff --git a/pyomo/contrib/edi/examples/aircraft_gp.py b/pyomo/contrib/edi/examples/aircraft_gp.py new file mode 100644 index 00000000000..f23da219b17 --- /dev/null +++ b/pyomo/contrib/edi/examples/aircraft_gp.py @@ -0,0 +1,123 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# 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. +# +# Development of this module was conducted as part of the Institute for +# the Design of Advanced Energy Systems (IDAES) with support through the +# Simulation-Based Engineering, Crosscutting Research Program within the +# U.S. Department of Energy’s Office of Fossil Energy and Carbon Management. +# +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +# =========== +# Description +# =========== +# A simple aircraft sizing problem, formulated as a Geometric Program +# From: Hoburg and Abbeel +# Geometric Programming for Aircraft Design Optimization +# AIAA Journal +# 2014 + +# ================= +# Import Statements +# ================= +import numpy as np +from pyomo.contrib.edi import Formulation + +# =================== +# Declare Formulation +# =================== +f = Formulation() + +# ================= +# Declare Variables +# ================= + +A = f.Variable(name="A", guess=10.0, units="-", description="aspect ratio") +C_D = f.Variable( + name="C_D", guess=0.025, units="-", description="Drag coefficient of wing" +) +C_f = f.Variable( + name="C_f", guess=0.003, units="-", description="skin friction coefficient" +) +C_L = f.Variable( + name="C_L", guess=0.5, units="-", description="Lift coefficient of wing" +) +D = f.Variable(name="D", guess=300, units="N", description="total drag force") +Re = f.Variable(name="Re", guess=3e6, units="-", description="Reynold's number") +S = f.Variable(name="S", guess=10.0, units="m^2", description="total wing area") +V = f.Variable(name="V", guess=30.0, units="m/s", description="cruising speed") +W = f.Variable(name="W", guess=10000.0, units="N", description="total aircraft weight") +W_w = f.Variable(name="W_w", guess=2500, units="N", description="wing weight") + +# ================= +# Declare Constants +# ================= +C_Lmax = f.Constant( + name="C_Lmax", value=2.0, units="-", description="max CL with flaps down" +) +CDA0 = f.Constant( + name="CDA0", value=0.0306, units="m^2", description="fuselage drag area" +) +e = f.Constant(name="e", value=0.96, units="-", description="Oswald efficiency factor") +k = f.Constant(name="k", value=1.2, units="-", description="form factor") +mu = f.Constant( + name="mu", value=1.78e-5, units="kg/m/s", description="viscosity of air" +) +N_ult = f.Constant( + name="N_ult", value=2.5, units="-", description="ultimate load factor" +) +rho = f.Constant(name="rho", value=1.23, units="kg/m^3", description="density of air") +S_wetratio = f.Constant( + name="Srat", value=2.05, units="-", description="wetted area ratio" +) +tau = f.Constant( + name="tau", value=0.12, units="-", description="airfoil thickness to chord ratio" +) +V_min = f.Constant(name="V_min", value=22, units="m/s", description="takeoff speed") +W_0 = f.Constant( + name="W_0", value=4940.0, units="N", description="aircraft weight excluding wing" +) +W_W_coeff1 = f.Constant( + name="W_c1", value=8.71e-5, units="1/m", description="Wing Weight Coefficient 1" +) +W_W_coeff2 = f.Constant( + name="W_c2", value=45.24, units="Pa", description="Wing Weight Coefficient 2" +) + +# ===================== +# Declare the Objective +# ===================== +f.Objective(D) + +# =================================== +# Declare some intermediate variables +# =================================== +pi = np.pi +C_D_fuse = CDA0 / S +C_D_wpar = k * C_f * S_wetratio +C_D_ind = C_L**2 / (pi * A * e) +W_w_strc = W_W_coeff1 * (N_ult * A**1.5 * (W_0 * W * S) ** 0.5) / tau +W_w_surf = W_W_coeff2 * S + +# ======================= +# Declare the Constraints +# ======================= +f.ConstraintList( + [ + C_D >= C_D_fuse + C_D_wpar + C_D_ind, + W_w >= W_w_surf + W_w_strc, + D >= 0.5 * rho * S * C_D * V**2, + Re == (rho / mu) * V * (S / A) ** 0.5, + C_f == 0.074 / Re**0.2, + W == 0.5 * rho * S * C_L * V**2, + W == 0.5 * rho * S * C_Lmax * V_min**2, + W >= W_0 + W_w, + ] +) diff --git a/pyomo/contrib/edi/examples/readme_example.py b/pyomo/contrib/edi/examples/readme_example.py new file mode 100644 index 00000000000..adf4799a569 --- /dev/null +++ b/pyomo/contrib/edi/examples/readme_example.py @@ -0,0 +1,97 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# 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. +# +# Development of this module was conducted as part of the Institute for +# the Design of Advanced Energy Systems (IDAES) with support through the +# Simulation-Based Engineering, Crosscutting Research Program within the +# U.S. Department of Energy’s Office of Fossil Energy and Carbon Management. +# +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +# ================= +# Import Statements +# ================= +import pyomo.environ as pyo +from pyomo.environ import units +from pyomo.contrib.edi import Formulation +from pyomo.contrib.edi import BlackBoxFunctionModel + +# =================== +# Declare Formulation +# =================== +f = Formulation() + +# ================= +# Declare Variables +# ================= +x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') +y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') +z = f.Variable(name='z', guess=1.0, units='m^2', description='Model output') + +# ================= +# Declare Constants +# ================= +c = f.Constant(name='c', value=1.0, units='', description='A constant c', size=2) + +# ===================== +# Declare the Objective +# ===================== +f.Objective(c[0] * x + c[1] * y) + + +# =================== +# Declare a Black Box +# =================== +class UnitCircle(BlackBoxFunctionModel): + def __init__(self): # The initialization function + # Initialize the black box model + super().__init__() + + # A brief description of the model + self.description = 'This model evaluates the function: z = x**2 + y**2' + + # Declare the black box model inputs + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + + # Declare the black box model outputs + self.outputs.append( + name='z', units='ft**2', description='Resultant of the unit circle' + ) + + # Declare the maximum available derivative + self.availableDerivative = 1 + + def BlackBox(self, x, y): # The actual function that does things + # Converts to correct units then casts to float + x = pyo.value(units.convert(x, self.inputs['x'].units)) + y = pyo.value(units.convert(y, self.inputs['y'].units)) + + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + + z *= self.outputs['z'].units + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + +# ======================= +# Declare the Constraints +# ======================= +f.ConstraintList([[z, '==', [x, y], UnitCircle()], z <= 1 * units.m**2]) + +# ================= +# Run the black box +# ================= +uc = UnitCircle() +bbo = uc.BlackBox(0.5 * units.m, 0.5 * units.m) diff --git a/pyomo/contrib/edi/formulation.py b/pyomo/contrib/edi/formulation.py new file mode 100644 index 00000000000..b1a2ae30f75 --- /dev/null +++ b/pyomo/contrib/edi/formulation.py @@ -0,0 +1,466 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# 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. +# +# Development of this module was conducted as part of the Institute for +# the Design of Advanced Energy Systems (IDAES) with support through the +# Simulation-Based Engineering, Crosscutting Research Program within the +# U.S. Department of Energy’s Office of Fossil Energy and Carbon Management. +# +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import pyomo +import pyomo.environ as pyo +from pyomo.util.check_units import assert_units_consistent +from pyomo.environ import ConcreteModel +from pyomo.environ import Var, Param, Objective, Constraint, Set +from pyomo.environ import maximize, minimize +from pyomo.environ import units as pyomo_units +from pyomo.common.dependencies import attempt_import + +egb, egb_available = attempt_import( + "pyomo.contrib.pynumero.interfaces.external_grey_box" +) +if egb_available: + from pyomo.contrib.pynumero.interfaces.external_grey_box import ( + ExternalGreyBoxModel, + ExternalGreyBoxBlock, + ) +else: + raise ImportError( + "pyomo.contrib.edi requires pyomo.contrib.pynumero to be installed to enable black box capability, this should have installed with base pyomo" + ) + +from pyomo.environ import ( + Reals, + PositiveReals, + NonPositiveReals, + NegativeReals, + NonNegativeReals, + Integers, + PositiveIntegers, + NonPositiveIntegers, + NegativeIntegers, + NonNegativeIntegers, + Boolean, + Binary, + Any, + AnyWithNone, + EmptySet, + UnitInterval, + PercentFraction, + RealInterval, + IntegerInterval, +) + +domainList = [ + Reals, + PositiveReals, + NonPositiveReals, + NegativeReals, + NonNegativeReals, + Integers, + PositiveIntegers, + NonPositiveIntegers, + NegativeIntegers, + NonNegativeIntegers, + Boolean, + Binary, + Any, + AnyWithNone, + EmptySet, + UnitInterval, + PercentFraction, + RealInterval, + IntegerInterval, +] + + +def decodeUnits(u_val): + if isinstance(u_val, str): + if u_val in ['', '-', 'None', ' ', 'dimensionless']: + return pyomo_units.__getattr__('dimensionless') + else: + return pyomo_units.__getattr__(u_val) + else: + return u_val + + +class Formulation(ConcreteModel): + def __init__(self): + super(Formulation, self).__init__() + # self._variable_counter = 1 + # self._constant_counter = 1 + self._objective_counter = 0 + self._constraint_counter = 0 + + self._variable_keys = [] + self._constant_keys = [] + self._objective_keys = [] + self._runtimeObjective_keys = [] + self._objective_keys = [] + self._runtimeConstraint_keys = [] + self._constraint_keys = [] + self._allConstraint_keys = [] + + def Variable( + self, name, guess, units, description='', size=None, bounds=None, domain=None + ): + if domain is None: + domain = Reals + else: + if domain not in domainList: + raise RuntimeError("Invalid domain") + + if bounds is not None: + if not isinstance(bounds, (list, tuple)): + raise ValueError( + 'The keyword bounds must be a 2 length list or tuple of floats' + ) + if len(bounds) != 2: + raise ValueError( + 'The keyword bounds must be a 2 length list or tuple of floats' + ) + if not isinstance(bounds[0], (float, int)): + raise ValueError( + 'The keyword bounds must be a 2 length list or tuple of floats' + ) + if not isinstance(bounds[1], (float, int)): + raise ValueError( + 'The keyword bounds must be a 2 length list or tuple of floats' + ) + if bounds[0] > bounds[1]: + raise ValueError("Lower bound is higher than upper bound") + + if size is not None: + if isinstance(size, (list, tuple)): + for i in range(0, len(size)): + if not isinstance(size[i], int): + raise ValueError( + 'Invalid size. Must be an integer or list/tuple of integers' + ) + if size[i] == 1 or size[i] == 0: + raise ValueError( + 'A value of 0 or 1 is not valid for defining size. Use fewer dimensions.' + ) + if i == 0: + st = pyo.Set(initialize=list(range(0, size[i]))) + else: + st *= pyo.Set(initialize=list(range(0, size[i]))) + st.construct() + self.add_component( + name, + pyo.Var( + st, + name=name, + initialize=guess, + domain=domain, + bounds=bounds, + doc=description, + units=decodeUnits(units), + ), + ) + else: + if isinstance(size, int): + if size == 1 or size == 0: + self.add_component( + name, + pyo.Var( + name=name, + initialize=guess, + domain=domain, + bounds=bounds, + doc=description, + units=decodeUnits(units), + ), + ) + else: + st = pyo.Set(initialize=list(range(0, size))) + st.construct() + self.add_component( + name, + pyo.Var( + st, + name=name, + initialize=guess, + domain=domain, + bounds=bounds, + doc=description, + units=decodeUnits(units), + ), + ) + else: + raise ValueError( + 'Invalid size. Must be an integer or list/tuple of integers' + ) + else: + self.add_component( + name, + pyo.Var( + name=name, + initialize=guess, + domain=domain, + bounds=bounds, + doc=description, + units=decodeUnits(units), + ), + ) + self.__dict__[name].construct() + self._variable_keys.append(name) + return self.__dict__[name] + + def Constant(self, name, value, units, description='', size=None, within=None): + if within is None: + within = Reals + else: + if within not in domainList: + raise RuntimeError("Invalid within") + + if size is not None: + if isinstance(size, (list, tuple)): + for i in range(0, len(size)): + if not isinstance(size[i], int): + raise ValueError( + 'Invalid size. Must be an integer or list/tuple of integers' + ) + if size[i] == 1: + raise ValueError( + 'A value of 1 is not valid for defining size. Use fewer dimensions.' + ) + if i == 0: + st = pyo.Set(initialize=list(range(0, size[i]))) + else: + st *= pyo.Set(initialize=list(range(0, size[i]))) + st.construct() + self.add_component( + name, + pyo.Param( + st, + name=name, + initialize=value, + within=within, + doc=description, + units=decodeUnits(units), + mutable=True, + ), + ) + else: + if isinstance(size, int): + if size == 1 or size == 0: + self.add_component( + name, + pyo.Param( + name=name, + initialize=value, + within=within, + doc=description, + units=decodeUnits(units), + mutable=True, + ), + ) + else: + st = pyo.Set(initialize=list(range(0, size))) + st.construct() + self.add_component( + name, + pyo.Param( + st, + name=name, + initialize=value, + within=within, + doc=description, + units=decodeUnits(units), + mutable=True, + ), + ) + else: + raise ValueError( + 'Invalid size. Must be an integer or list/tuple of integers' + ) + else: + self.add_component( + name, + pyo.Param( + name=name, + initialize=value, + within=within, + doc=description, + units=decodeUnits(units), + mutable=True, + ), + ) + + self.__dict__[name].construct() + self._constant_keys.append(name) + return self.__dict__[name] + + def Objective(self, expr, sense=minimize): + self._objective_counter += 1 + self.add_component( + 'objective_' + str(self._objective_counter), + pyo.Objective(expr=expr, sense=sense), + ) + self._objective_keys.append('objective_' + str(self._objective_counter)) + self.__dict__['objective_' + str(self._objective_counter)].construct() + + # def RuntimeObjective(self): + # pass + + def Constraint(self, expr): + self._constraint_counter += 1 + conName = 'constraint_' + str(self._constraint_counter) + self.add_component(conName, pyo.Constraint(expr=expr)) + self._constraint_keys.append(conName) + self._allConstraint_keys.append(conName) + self.__dict__[conName].construct() + + def RuntimeConstraint(self, outputs, operators, inputs, black_box): + self._constraint_counter += 1 + conName = 'constraint_' + str(self._constraint_counter) + self._runtimeConstraint_keys.append(conName) + self._allConstraint_keys.append(conName) + + self.add_component(conName, ExternalGreyBoxBlock()) + self.__dict__[conName].construct() + + # TODO: Need to include operators after Michael fixes things + + inputs_raw = inputs + outputs_raw = outputs + operators_raw = operators + + if isinstance( + inputs_raw, (pyomo.core.base.var.IndexedVar, pyomo.core.base.var.ScalarVar) + ): + inputs_raw = [inputs_raw] + elif isinstance(inputs_raw, (list, tuple)): + inputs_raw = list(inputs_raw) + else: + raise ValueError("Invalid type for input variables") + + if isinstance( + outputs_raw, (pyomo.core.base.var.IndexedVar, pyomo.core.base.var.ScalarVar) + ): + outputs_raw = [outputs_raw] + elif isinstance(outputs_raw, (list, tuple)): + outputs_raw = list(outputs_raw) + else: + raise ValueError("Invalid type for output variables") + for lst in [outputs_raw, inputs_raw]: + for vr in lst: + if not isinstance( + vr, (pyomo.core.base.var.IndexedVar, pyomo.core.base.var.ScalarVar) + ): + raise ValueError("Invalid type when checking inputs and outputs") + + if isinstance(operators_raw, (list, tuple)): + operators_raw = list(operators_raw) + elif isinstance(operators_raw, str): + operators_raw = [operators_raw] + else: + raise ValueError("Invalid type for operators") + for opr in operators_raw: + if opr not in ["==", ">=", "<="]: + raise ValueError("Invalid operator") + + black_box.setOptimizationVariables(inputs_raw, outputs_raw) + + outputs_raw_length = len(outputs_raw) + operators_raw_length = len(operators_raw) + + outputs_unwrapped = [] + for ovar in outputs_raw: + if isinstance(ovar, pyomo.core.base.var.ScalarVar): + outputs_unwrapped.append(ovar) + else: # isinstance(ovar, pyomo.core.base.var.IndexedVar), validated above + validIndices = list(ovar.index_set().data()) + for vi in validIndices: + outputs_unwrapped.append(ovar[vi]) + + inputs_unwrapped = [] + for ivar in inputs_raw: + if isinstance(ivar, pyomo.core.base.var.ScalarVar): + inputs_unwrapped.append(ivar) + else: # isinstance(ivar, pyomo.core.base.var.IndexedVar), validated above + validIndices = list(ivar.index_set().data()) + for vi in validIndices: + inputs_unwrapped.append(ivar[vi]) + + black_box._NunwrappedOutputs = len(outputs_unwrapped) + black_box._NunwrappedInputs = len(inputs_unwrapped) + black_box.post_init_setup() + + # TODO: Need to unwrap operators + + self.__dict__[conName].set_external_model( + black_box, inputs=inputs_unwrapped, outputs=outputs_unwrapped + ) # ,operators=operators_unwrapped) + + def ConstraintList(self, conList): + for i in range(0, len(conList)): + con = conList[i] + if isinstance(con, (tuple, list)): + self.RuntimeConstraint(*con) + elif isinstance(con, dict): + self.RuntimeConstraint(**con) + else: + self.Constraint(con) + + def get_variables(self): + return [ + self.__dict__[nm] + for nm in self.__dict__.keys() + if nm in self._variable_keys + ] + + def get_constants(self): + return [ + self.__dict__[nm] + for nm in self.__dict__.keys() + if nm in self._constant_keys + ] + + def get_objectives(self): + return [ + self.__dict__[nm] + for nm in self.__dict__.keys() + if nm in self._objective_keys + ] + + def get_constraints(self): + return [ + self.__dict__[nm] + for nm in self.__dict__.keys() + if nm in self._allConstraint_keys + ] + + def get_explicitConstraints(self): + return [ + self.__dict__[nm] + for nm in self.__dict__.keys() + if nm in self._constraint_keys + ] + + def get_runtimeConstraints(self): + return [ + self.__dict__[nm] + for nm in self.__dict__.keys() + if nm in self._runtimeConstraint_keys + ] + + def check_units(self): + for i in range(1, self._objective_counter + 1): + assert_units_consistent(self.__dict__['objective_' + str(i)]) + + for i in range(1, self._constraint_counter + 1): + if not isinstance( + self.__dict__['constraint_' + str(i)], + pyomo.contrib.pynumero.interfaces.external_grey_box.ExternalGreyBoxBlock, + ): + assert_units_consistent(self.__dict__['constraint_' + str(i)]) diff --git a/pyomo/contrib/edi/tests/__init__.py b/pyomo/contrib/edi/tests/__init__.py new file mode 100644 index 00000000000..38b30839be3 --- /dev/null +++ b/pyomo/contrib/edi/tests/__init__.py @@ -0,0 +1,16 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# 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. +# +# Development of this module was conducted as part of the Institute for +# the Design of Advanced Energy Systems (IDAES) with support through the +# Simulation-Based Engineering, Crosscutting Research Program within the +# U.S. Department of Energy’s Office of Fossil Energy and Carbon Management. +# +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ diff --git a/pyomo/contrib/edi/tests/test_blackbox.py b/pyomo/contrib/edi/tests/test_blackbox.py new file mode 100644 index 00000000000..7f273763317 --- /dev/null +++ b/pyomo/contrib/edi/tests/test_blackbox.py @@ -0,0 +1,1725 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# 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. +# +# Development of this module was conducted as part of the Institute for +# the Design of Advanced Energy Systems (IDAES) with support through the +# Simulation-Based Engineering, Crosscutting Research Program within the +# U.S. Department of Energy’s Office of Fossil Energy and Carbon Management. +# +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import pyomo.common.unittest as unittest +import pyomo.environ as pyo +from pyomo.common.dependencies import attempt_import + +from pyomo.core.base.units_container import pint_available + +from pyomo.common.dependencies import numpy, numpy_available +from pyomo.common.dependencies import scipy, scipy_available + +egb, egb_available = attempt_import( + "pyomo.contrib.pynumero.interfaces.external_grey_box" +) + +formulation_available = False +try: + from pyomo.contrib.edi import Formulation + + formulation_available = True +except: + pass + # formulation_available = False + +blackbox_available = False +try: + from pyomo.contrib.edi import BlackBoxFunctionModel + + blackbox_available = True +except: + pass + # blackbox_available = False + + +if numpy_available: + import numpy as np + + +@unittest.skipIf( + not egb_available, 'Testing pyomo.contrib.edi requires pynumero external grey boxes' +) +@unittest.skipIf(not formulation_available, 'Formulation import failed') +@unittest.skipIf(not blackbox_available, 'Blackbox import failed') +@unittest.skipIf(not numpy_available, 'Testing pyomo.contrib.edi requires numpy') +@unittest.skipIf(not scipy_available, 'Testing pyomo.contrib.edi requires scipy') +@unittest.skipIf(not pint_available, 'Testing units requires pint') +class TestEDIBlackBox(unittest.TestCase): + def test_edi_blackbox_variable(self): + "Tests the black box variable class" + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + x = BlackBoxFunctionModel_Variable('x', '') + x_print = x.__repr__() + x_name = x.name + x_units = x.units + x_size = x.size + x_desc = x.description + self.assertRaises(ValueError, x.__init__, *(1.0, '')) + + x.__init__('x', units.dimensionless) + + self.assertRaises(ValueError, x.__init__, *('x', 1.0)) + + x.__init__('x', units.dimensionless, '', 'flex') + x.__init__('x', units.dimensionless, '', ['flex', 2]) + x.__init__('x', units.dimensionless, '', 2) + x.__init__('x', units.dimensionless, '', [2, 2]) + + self.assertRaises(ValueError, x.__init__, *('x', '', '', [[], 2])) + self.assertRaises(ValueError, x.__init__, *('x', '', '', [2, 1])) + + x.__init__('x', units.dimensionless, '', None) + x.__init__('x', None, '', None) + + self.assertRaises(ValueError, x.__init__, *('x', '', '', 1)) + self.assertRaises(ValueError, x.__init__, *('x', '', '', {})) + self.assertRaises(ValueError, x.__init__, *('x', '', 1.0)) + + def test_edi_blackbox_tcl(self): + "Tests the black box type checked list class" + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + tcl = TypeCheckedList(int, [1, 2, 3]) + tcl[1] = 1 + tcl[0:2] = [1, 2] + + self.assertRaises(ValueError, tcl.__init__, *(int, 1)) + self.assertRaises(ValueError, tcl.__setitem__, *(1, 3.333)) + self.assertRaises(ValueError, tcl.__setitem__, *(1, [1, 2.222])) + self.assertRaises(ValueError, tcl.append, *(2.222,)) + + def test_edi_blackbox_bbl(self): + "Tests the black box BBList class" + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + bbl = BBList() + bbl.append('x', '') + bbl.append('y', '') + bbl.append(BlackBoxFunctionModel_Variable('z', '')) + bbl.append(var=BlackBoxFunctionModel_Variable('u', '')) + self.assertRaises( + ValueError, bbl.append, *(BlackBoxFunctionModel_Variable('u', ''),) + ) + self.assertRaises(ValueError, bbl.append, *('badvar',)) + self.assertRaises(ValueError, bbl.append, *(2.222,)) + + self.assertRaises(ValueError, bbl.append, *('bv', '', ''), **{'units': 'm'}) + self.assertRaises(ValueError, bbl.append, **{'units': 'm', 'description': 'hi'}) + self.assertRaises(ValueError, bbl.append, **{'name': 'x', 'units': ''}) + self.assertRaises(ValueError, bbl.append, *('bv', '', '', 0, 'extra')) + + xv = bbl['x'] + xv2 = bbl[0] + self.assertRaises(ValueError, bbl.__getitem__, *(2.22,)) + + def test_edi_blackbox_someexceptions(self): + "Tests some of the exceptions in the black box model class" + import numpy as np + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + bb = BlackBoxFunctionModel() + bb.inputVariables_optimization = [1, 2, 3] + # bb.set_input_values(np.array([1,2,3])) + self.assertRaises(ValueError, bb.input_names, *()) + # self.assertRaises(ValueError,bb.fillCache,*( )) + + bb = BlackBoxFunctionModel() + bb.outputVariables_optimization = [1, 2, 3] + self.assertRaises(ValueError, bb.output_names, *()) + + def test_edi_blackbox_etc_1(self): + "Tests a black box assertion issue" + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + bbfm = BlackBoxFunctionModel() + self.assertRaises(AttributeError, bbfm.BlackBox, ()) + + def test_edi_blackbox_etc_2(self): + "Tests a black box assertion issue" + import numpy as np + from pyomo.environ import units + import pyomo.environ as pyo + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='y variable') + z = f.Variable(name='z', guess=1.0, units='m^2', description='Output var') + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', units='ft**2', description='Output variable' + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + # Converts to correct units then casts to float + x = pyo.value(units.convert(x, self.inputs[0].units)) + y = pyo.value(units.convert(y, self.inputs[1].units)) + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + f.ConstraintList([z <= 1 * units.m**2, [z, '==', [x, y], UnitCircle()]]) + + f.__dict__['constraint_2'].get_external_model().set_input_values( + np.array([2.0, 2.0]) + ) + f.__dict__['constraint_2'].get_external_model().inputVariables_optimization = [ + 1, + 2, + ] + # f.__dict__['constraint_2'].get_external_model().fillCache() + self.assertRaises( + ValueError, f.__dict__['constraint_2'].get_external_model().fillCache, *() + ) + + def test_edi_blackbox_etc_3(self): + "Tests a black box assertion issue" + import numpy as np + from pyomo.environ import units + import pyomo.environ as pyo + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='y variable') + z = f.Variable(name='z', guess=1.0, units='m^2', description='Output var') + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', units='ft**2', description='Output variable' + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + # Converts to correct units then casts to float + x = pyo.value(units.convert(x, self.inputs[0].units)) + y = pyo.value(units.convert(y, self.inputs[1].units)) + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + return ['err'], [[dzdx, dzdy]] # return z, grad(z), hess(z)... + + f.ConstraintList([z <= 1 * units.m**2, [z, '==', [x, y], UnitCircle()]]) + + self.assertRaises( + ValueError, f.__dict__['constraint_2'].get_external_model().fillCache, *() + ) + + def test_edi_blackbox_etc_4(self): + "Tests a black box assertion issue" + import numpy as np + from pyomo.environ import units + import pyomo.environ as pyo + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='y variable') + z = f.Variable(name='z', guess=1.0, units='m^2', description='Output var') + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', units='ft**2', description='Output variable' + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + # Converts to correct units then casts to float + x = pyo.value(units.convert(x, self.inputs[0].units)) + y = pyo.value(units.convert(y, self.inputs[1].units)) + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + f.ConstraintList([z <= 1 * units.m**2, [z, '==', [x, y], UnitCircle()]]) + + f.__dict__['constraint_2'].get_external_model().set_input_values( + np.array([2.0, 2.0]) + ) + f.__dict__['constraint_2'].get_external_model().outputVariables_optimization = [ + 1 + ] + # f.__dict__['constraint_2'].get_external_model().fillCache() + self.assertRaises( + ValueError, f.__dict__['constraint_2'].get_external_model().fillCache, *() + ) + + def test_edi_blackbox_example_1(self): + "Tests a black box example construction" + from pyomo.environ import units + import pyomo.environ as pyo + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='y variable') + z = f.Variable(name='z', guess=1.0, units='m^2', description='Output var') + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', units='ft**2', description='Output variable' + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + # Converts to correct units then casts to float + x = pyo.value(units.convert(x, self.inputs[0].units)) + y = pyo.value(units.convert(y, self.inputs[1].units)) + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + f.ConstraintList([z <= 1 * units.m**2, [z, '==', [x, y], UnitCircle()]]) + + f.__dict__['constraint_2'].get_external_model().set_input_values( + np.array([2.0, 2.0]) + ) + opt = f.__dict__['constraint_2'].get_external_model().evaluate_outputs() + jac = ( + f.__dict__['constraint_2'] + .get_external_model() + .evaluate_jacobian_outputs() + .todense() + ) + + self.assertAlmostEqual(opt[0], 8) + self.assertAlmostEqual(jac[0, 0], 4) + self.assertAlmostEqual(jac[0, 1], 4) + + sm = f.__dict__['constraint_2'].get_external_model().summary + e_print = f.__dict__['constraint_2'].get_external_model().__repr__() + + def test_edi_blackbox_example_2(self): + "Tests a black box example construction" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable( + name='x', guess=1.0, units='m', description='The x variable', size=3 + ) + y = f.Variable( + name='y', guess=1.0, units='m**2', description='The y variable', size=3 + ) + f.Objective(y[0] + y[1] + y[2]) + + class Parabola(BlackBoxFunctionModel): + def __init__(self): + super(Parabola, self).__init__() + self.description = 'This model evaluates the function: y = x**2' + self.inputs.append( + name='x', size=3, units='ft', description='The x variable' + ) + self.outputs.append( + name='y', size=3, units='ft**2', description='The y variable' + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, *args, **kwargs): # The actual function that does things + runCases, remainingKwargs = self.parseInputs(*args, **kwargs) + + x = self.sanitizeInputs(runCases[0]['x']) + x = np.array([pyo.value(xval) for xval in x], dtype=np.float64) + + y = x**2 # Compute y + dydx = 2 * x # Compute dy/dx + + y = y * units.ft**2 + dydx = np.diag(dydx) + dydx = dydx * units.ft # units.ft**2 / units.ft + + return y, [dydx] # return z, grad(z), hess(z)... + + f.ConstraintList( + [{'outputs': y, 'operators': '==', 'inputs': x, 'black_box': Parabola()}] + ) + + f.__dict__['constraint_1'].get_external_model().set_input_values(np.ones(3) * 2) + opt = f.__dict__['constraint_1'].get_external_model().evaluate_outputs() + jac = ( + f.__dict__['constraint_1'] + .get_external_model() + .evaluate_jacobian_outputs() + .todense() + ) + + self.assertAlmostEqual(opt[0], 4) + self.assertAlmostEqual(opt[1], 4) + self.assertAlmostEqual(opt[2], 4) + self.assertAlmostEqual(jac[0, 0], 4) + self.assertAlmostEqual(jac[0, 1], 0) + self.assertAlmostEqual(jac[0, 2], 0) + self.assertAlmostEqual(jac[0, 1], 0) + self.assertAlmostEqual(jac[1, 1], 4) + self.assertAlmostEqual(jac[2, 1], 0) + self.assertAlmostEqual(jac[2, 0], 0) + self.assertAlmostEqual(jac[2, 1], 0) + self.assertAlmostEqual(jac[2, 2], 4) + + sm = f.__dict__['constraint_1'].get_external_model().summary + e_print = f.__dict__['constraint_1'].get_external_model().__repr__() + + def test_edi_blackbox_example_3(self): + "Tests a black box example construction" + import numpy as np + from pyomo.environ import units + import pyomo.environ as pyo + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='', description='x variable', size=3) + y = f.Variable(name='y', guess=1.0, units='', description='y variable') + f.Objective(y) + + class Norm_2(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the two norm' + self.inputs.append( + name='x', units='', description='The x variable', size=3 + ) + self.outputs.append(name='y', units='', description='The y variable') + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, *args, **kwargs): # The actual function that does things + runCases, remainingKwargs = self.parseInputs(*args, **kwargs) + + x = self.sanitizeInputs(runCases[0]['x']) + x = np.array([pyo.value(xval) for xval in x], dtype=np.float64) + + y = x[0] ** 2 + x[1] ** 2 + x[2] ** 2 # Compute y + dydx0 = 2 * x[0] # Compute dy/dx0 + dydx1 = 2 * x[1] # Compute dy/dx1 + dydx2 = 2 * x[2] # Compute dy/dx2 + + y = y * units.dimensionless + dydx = np.array([dydx0, dydx1, dydx2]) * units.dimensionless + # dydx0 = dydx0 * units.dimensionless + # dydx1 = dydx1 * units.dimensionless + # dydx2 = dydx2 * units.dimensionless + + return y, [dydx] # return z, grad(z), hess(z)... + + f.ConstraintList( + [{'outputs': y, 'operators': '==', 'inputs': x, 'black_box': Norm_2()}] + ) + + f.__dict__['constraint_1'].get_external_model().set_input_values(np.ones(3) * 2) + opt = f.__dict__['constraint_1'].get_external_model().evaluate_outputs() + jac = ( + f.__dict__['constraint_1'] + .get_external_model() + .evaluate_jacobian_outputs() + .todense() + ) + + self.assertAlmostEqual(opt[0], 12) + self.assertAlmostEqual(jac[0, 0], 4) + self.assertAlmostEqual(jac[0, 1], 4) + self.assertAlmostEqual(jac[0, 2], 4) + + sm = f.__dict__['constraint_1'].get_external_model().summary + e_print = f.__dict__['constraint_1'].get_external_model().__repr__() + + def test_edi_blackbox_example_4(self): + "Tests a black box example construction" + import numpy as np + from pyomo.environ import units + import pyomo.environ as pyo + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='', description='x variable') + y = f.Variable(name='y', guess=1.0, units='', description='y variable', size=3) + f.Objective(y[0] ** 2 + y[1] ** 2 + y[2] ** 2) + + class VectorCast(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the two norm' + self.inputs.append(name='x', units='', description='The x variable') + self.outputs.append( + name='y', units='', description='The y variable', size=3 + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, *args, **kwargs): # The actual function that does things + runCases, remainingKwargs = self.parseInputs(*args, **kwargs) + + x = pyo.value(self.sanitizeInputs(runCases[0]['x'])) + + y = np.array([x, x, x]) * units.dimensionless + dydx = np.array([1.0, 1.0, 1.0]) * units.dimensionless + + return y, [dydx] # return z, grad(z), hess(z)... + + f.ConstraintList( + [{'outputs': y, 'operators': '==', 'inputs': x, 'black_box': VectorCast()}] + ) + + f.__dict__['constraint_1'].get_external_model().set_input_values(np.ones(3) * 2) + opt = f.__dict__['constraint_1'].get_external_model().evaluate_outputs() + jac = ( + f.__dict__['constraint_1'] + .get_external_model() + .evaluate_jacobian_outputs() + .todense() + ) + + self.assertAlmostEqual(opt[0], 2.0) + self.assertAlmostEqual(opt[1], 2.0) + self.assertAlmostEqual(opt[2], 2.0) + self.assertAlmostEqual(jac[0, 0], 1.0) + self.assertAlmostEqual(jac[1, 0], 1.0) + self.assertAlmostEqual(jac[2, 0], 1.0) + + sm = f.__dict__['constraint_1'].get_external_model().summary + e_print = f.__dict__['constraint_1'].get_external_model().__repr__() + + def test_edi_blackbox_badexample_1(self): + "Tests a black box example construction" + import numpy as np + from pyomo.environ import units + import pyomo.environ as pyo + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='', description='x variable', size=3) + y = f.Variable(name='y', guess=1.0, units='', description='y variable') + f.Objective(y) + + class Norm_2(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the two norm' + self.inputs.append( + name='x', units='', description='The x variable', size=3 + ) + self.outputs.append(name='y', units='', description='The y variable') + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, *args, **kwargs): # The actual function that does things + runCases, remainingKwargs = self.parseInputs(*args, **kwargs) + + x = self.sanitizeInputs(runCases[0]['x']) + x = np.array([pyo.value(xval) for xval in x], dtype=np.float64) + + y = x[0] ** 2 + x[1] ** 2 + x[2] ** 2 # Compute y + dydx0 = 2 * x[0] # Compute dy/dx0 + dydx1 = 2 * x[1] # Compute dy/dx1 + dydx2 = 2 * x[2] # Compute dy/dx2 + + y = y * units.dimensionless + # dydx = np.array([dydx0,dydx1,dydx2]) * units.dimensionless + dydx0 = dydx0 * units.dimensionless + dydx1 = dydx1 * units.dimensionless + dydx2 = dydx2 * units.dimensionless + + return y, [[dydx0, dydx1, dydx2]] # return z, grad(z), hess(z)... + + f.ConstraintList( + [{'outputs': y, 'operators': '==', 'inputs': x, 'black_box': Norm_2()}] + ) + + f.__dict__['constraint_1'].get_external_model().set_input_values(np.ones(3) * 2) + self.assertRaises( + ValueError, + f.__dict__['constraint_1'].get_external_model().evaluate_outputs, + *() + ) + + def test_edi_blackbox_smallfunctions(self): + "Tests the more general value and convert functions" + import numpy as np + from pyomo.environ import units + import pyomo.environ as pyo + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + bb = BlackBoxFunctionModel() + t1 = bb.convert(2 * units.m, units.ft) + t2 = bb.convert(np.ones([2, 2]) * units.m, units.ft) + self.assertRaises(ValueError, bb.convert, *('err', units.ft)) + + t3 = bb.pyomo_value(2 * units.m) + t3 = bb.pyomo_value(np.ones([2, 2]) * units.m) + t4 = bb.pyomo_value(np.ones([2, 2])) + + bb.sizeCheck([2, 2], np.ones([2, 2]) * units.m) + self.assertRaises(ValueError, bb.sizeCheck, *(2, 3 * units.m)) + self.assertRaises( + ValueError, bb.sizeCheck, *([10, 10, 10], np.ones([2, 2]) * units.m) + ) + self.assertRaises( + ValueError, bb.sizeCheck, *([10, 10], np.ones([2, 2]) * units.m) + ) + self.assertRaises(ValueError, bb.sizeCheck, *([10, 10], [])) + + def test_edi_blackbox_bare_example_1(self): + "Tests a black box example construction without an optimization problem" + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + class SignomialTest(BlackBoxFunctionModel): + def __init__(self): + # Set up all the attributes by calling Model.__init__ + super().__init__() + + # Setup Inputs + self.inputs.append('x', '', 'Independent Variable') + + # Setup Outputs + self.outputs.append('y', '', 'Dependent Variable') + + # Simple model description + self.description = ( + 'This model evaluates the function: max([-6*x-6, x**4-3*x**2])' + ) + + self.availableDerivative = 1 + + # standard function call is y(, dydx, ...) = self.BlackBox(**{'x1':x1, 'x2':x2, ...}) + def BlackBox(self, *args, **kwargs): + runCases, extras = self.parseInputs(*args, **kwargs) + x = np.array( + [pyo.value(runCases[i]['x']) for i in range(0, len(runCases))] + ) + + y = np.maximum(-6 * x - 6, x**4 - 3 * x**2) + dydx = 4 * x**3 - 6 * x + ddy_ddx = 12 * x**2 - 6 + gradientSwitch = -6 * x - 6 > x**4 - 3 * x**2 + dydx[gradientSwitch] = -6 + ddy_ddx[gradientSwitch] = 0 + + y = [yval * units.dimensionless for yval in y] + dydx = [dydx[i] * units.dimensionless for i in range(0, len(dydx))] + + returnMode = 1 + + if returnMode < 0: + returnMode = -1 * (returnMode + 1) + if returnMode == 0: + return y[0] + if returnMode == 1: + return y[0], dydx + else: + if returnMode == 0: + opt = [] + for i in range(0, len(y)): + opt.append([y[i]]) + return opt + if returnMode == 1: + opt = [] + for i in range(0, len(y)): + opt.append([[y[i]], [[[dydx[i]]]]]) + return opt + + s = SignomialTest() + ivals = [[x] for x in np.linspace(-2, 2, 11)] + + # How the black box may be called using EDI + bbo = s.BlackBox(**{'x': 0.5}) + bbo = s.BlackBox({'x': 0.5}) + bbo = s.BlackBox(**{'x': 0.5, 'optn': True}) + + # Additional options available with parseInputs + bbo = s.BlackBox(*[0.5], **{'optn1': True, 'optn2': False}) + bbo = s.BlackBox(*[0.5, True], **{'optn': False}) + bbo = s.BlackBox({'x': [x for x in np.linspace(-2, 2, 11)]}) + bbo = s.BlackBox([{'x': x} for x in np.linspace(-2, 2, 11)]) + bbo = s.BlackBox([[x] for x in np.linspace(-2, 2, 11)]) + bbo = s.BlackBox([[x] for x in np.linspace(-2, 2, 11)], True, optn=False) + bbo = s.BlackBox([[x] for x in np.linspace(-2, 2, 11)], optn1=True, optn2=False) + + sm = s.summary + + self.assertRaises(ValueError, s.BlackBox, *([[[1, 2], 2, 3], [1, 2, 3]])) + self.assertRaises(ValueError, s.BlackBox, *([[np.ones(2), 2, 3], [1, 2, 3]])) + + self.assertRaises(ValueError, s.sanitizeInputs, *()) + self.assertRaises(ValueError, s.sanitizeInputs, *(1, 2, 3)) + self.assertRaises( + ValueError, s.sanitizeInputs, **{'z': 2.0 * units.dimensionless} + ) + self.assertRaises(ValueError, s.sanitizeInputs, *(2.0 * units.ft,)) + self.assertRaises(NotImplementedError, s.checkOutputs, *()) + + def test_edi_blackbox_bare_example_2(self): + "Tests a black box example construction without an optimization problem" + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + from pyomo.common.formatting import tostr + + class PassThrough(BlackBoxFunctionModel): + def __init__(self): + # Set up all the attributes by calling Model.__init__ + super().__init__() + + # Setup Inputs + self.inputs.append('x', '', 'Independent Variable', size=[2, 2]) + + # Setup Outputs + self.outputs.append('y', '', 'Dependent Variable', size=[2, 2]) + + # Simple model description + self.description = 'This model is a pass through)' + + self.availableDerivative = 1 + + # standard function call is y(, dydx, ...) = self.BlackBox(**{'x1':x1, 'x2':x2, ...}) + def BlackBox(self, *args, **kwargs): + runCases, extras = self.parseInputs(*args, **kwargs) + x = [ + self.pyomo_value(runCases[i]['x']) for i in range(0, len(runCases)) + ] + + y = [] + dydx = [] + + for xval in x: + y.append(xval * units.dimensionless) + dydx_temp = np.zeros([2, 2, 2, 2]) + dydx_temp[0, 0, 0, 0] = 1.0 + dydx_temp[0, 1, 0, 1] = 1.0 + dydx_temp[1, 0, 1, 0] = 1.0 + dydx_temp[1, 1, 1, 1] = 1.0 + + dydx.append(dydx_temp * units.dimensionless) + + returnMode = 1 + + if returnMode < 0: + returnMode = -1 * (returnMode + 1) + if returnMode == 0: + return y[0] + if returnMode == 1: + return y[0], dydx + else: + if returnMode == 0: + opt = [] + for i in range(0, len(y)): + opt.append([y[i]]) + return opt + if returnMode == 1: + opt = [] + for i in range(0, len(y)): + opt.append([[y[i]], [[[dydx[i]]]]]) + return opt + + bb = PassThrough() + ivals = [ + [np.eye(2) * units.dimensionless], + [np.ones([2, 2]) * units.dimensionless], + [np.zeros([2, 2]) * units.dimensionless], + ] + + xv = np.eye(2) * units.dimensionless + + # How the black box may be called using EDI + bbo = bb.BlackBox(**{'x': xv}) + bbo = bb.BlackBox({'x': xv}) + bbo = bb.BlackBox(**{'x': xv, 'optn': True}) + + # # Additional options available with parseInputs + bbo = bb.BlackBox(*[xv], **{'optn1': True, 'optn2': False}) + bbo = bb.BlackBox(*[xv, True], **{'optn': False}) + bbo = bb.BlackBox({'x': [x[0] for x in ivals]}) + bbo = bb.BlackBox([{'x': x[0]} for x in ivals]) + bbo = bb.BlackBox([[x[0]] for x in ivals]) + bbo = bb.BlackBox([[x[0]] for x in ivals], True, optn=False) + bbo = bb.BlackBox([[x[0]] for x in ivals], optn1=True, optn2=False) + + sm = bb.summary + + self.assertRaises(ValueError, bb.sanitizeInputs, *(np.ones([2, 2]) * units.ft,)) + + def test_edi_blackbox_bare_example_3(self): + "Tests a black box example construction" + import numpy as np + from pyomo.environ import units + import pyomo.environ as pyo + from pyomo.contrib.edi import Formulation + from pyomo.contrib.edi.blackBoxFunctionModel import ( + BlackBoxFunctionModel_Variable, + TypeCheckedList, + BBList, + BlackBoxFunctionModel, + ) + + class Norm_2(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the two norm' + self.inputs.append( + name='x', units='', description='The x variable', size=3 + ) + self.inputs.append( + name='y', units='', description='The y variable', size=2 + ) + self.outputs.append(name='z', units='', description='The z variable') + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, *args, **kwargs): # The actual function that does things + runCases, remainingKwargs = self.parseInputs(*args, **kwargs) + + x = [rc['x'] for rc in runCases][0] + x = np.array([self.pyomo_value(xval) for xval in x], dtype=np.float64) + + y = [rc['y'] for rc in runCases][0] + y = np.array([self.pyomo_value(yval) for yval in y], dtype=np.float64) + + z = x[0] ** 2 + x[1] ** 2 + x[2] ** 2 + y[0] ** 2 + y[1] ** 2 + dzdx0 = 2 * x[0] # Compute dy/dx0 + dzdx1 = 2 * x[1] # Compute dy/dx1 + dzdx2 = 2 * x[2] # Compute dy/dx2 + dzdy0 = 2 * y[0] + dzdy1 = 2 * y[1] + + z = z * units.dimensionless + dz = np.array([dzdx0, dzdx1, dzdx2, dzdy0, dzdy1]) * units.dimensionless + # dydx0 = dydx0 * units.dimensionless + # dydx1 = dydx1 * units.dimensionless + # dydx2 = dydx2 * units.dimensionless + + return z, [dz] # return z, grad(z), hess(z)... + + bb = Norm_2() + bbo = bb.BlackBox( + { + 'x': np.array([0, 0, 0]) * units.dimensionless, + 'y': np.array([0, 0]) * units.dimensionless, + } + ) + bbo = bb.BlackBox( + np.array([0, 0, 0]) * units.dimensionless, + y=np.array([0, 0]) * units.dimensionless, + ) + + self.assertRaises( + ValueError, + bb.BlackBox, + *( + { + 'er': np.array([0, 0, 0]) * units.dimensionless, + 'y': np.array([0, 0]) * units.dimensionless, + }, + ) + ) + self.assertRaises(ValueError, bb.BlackBox, *('err',)) + self.assertRaises( + ValueError, + bb.BlackBox, + *( + np.array([0, 0, 0]) * units.dimensionless, + np.array([0, 0]) * units.dimensionless, + ), + **{'x': 'err too many'} + ) + self.assertRaises( + ValueError, + bb.BlackBox, + *(np.array([0, 0, 0]) * units.dimensionless,), + **{'notY': np.array([0, 0]) * units.dimensionless} + ) + + def test_edi_blackbox_bare_example_4(self): + "Tests a black box example construction without an optimization problem" + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + from pyomo.common.formatting import tostr + + class PassThrough(BlackBoxFunctionModel): + def __init__(self): + # Set up all the attributes by calling Model.__init__ + super().__init__() + + # Setup Inputs + self.inputs.append('x', '', 'X Variable') + self.inputs.append('y', '', 'Y Variable') + + # Setup Outputs + self.outputs.append('u', '', 'U Variable') + self.outputs.append('v', '', 'V Variable') + + # Simple model description + self.description = 'This model is a pass through)' + + self.availableDerivative = 1 + + # standard function call is y(, dydx, ...) = self.BlackBox(**{'x1':x1, 'x2':x2, ...}) + def BlackBox(self, *args, **kwargs): + runCases, extras = self.parseInputs(*args, **kwargs) + x = [ + self.pyomo_value(runCases[i]['x']) for i in range(0, len(runCases)) + ] + y = [ + self.pyomo_value(runCases[i]['y']) for i in range(0, len(runCases)) + ] + + u = [] + dudx = [] + dudy = [] + v = [] + dvdx = [] + dvdy = [] + + for xval in x: + u.append(xval * units.dimensionless) + dudx.append(1.0 * units.dimensionless) + dudy.append(0.0 * units.dimensionless) + + for yval in y: + v.append(yval * units.dimensionless) + dvdx.append(0.0 * units.dimensionless) + dvdy.append(1.0 * units.dimensionless) + + returnMode = 1 + + if returnMode < 0: + returnMode = -1 * (returnMode + 1) + if returnMode == 0: + return [u[0], v[0]] + if returnMode == 1: + return [u[0], v[0]], [[dudx[0], dudy[0]], [dvdx[0], dvdy[0]]] + else: + if returnMode == 0: + opt = [] + for i in range(0, len(y)): + opt.append([u[i], v[i]]) + return opt + if returnMode == 1: + opt = [] + for i in range(0, len(y)): + opt.append( + [[u[i], v[i]], [[dudx[i], dudy[i]], [dvdx[i], dvdy[i]]]] + ) + return opt + + bb = PassThrough() + bbo = bb.BlackBox(1.0, 1.0) + bbo = bb.BlackBox({'x': np.linspace(0, 10, 11), 'y': np.linspace(0, 10, 11)}) + + def test_edi_blackbox_standardized_1(self): + "Tests the blackbox_standardized function" + import pyomo + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel, ediprint + from pyomo.common.formatting import tostr + from pyomo.environ import units as pyomo_units + + class PassThrough(BlackBoxFunctionModel): + def __init__(self): + # Set up all the attributes by calling Model.__init__ + super().__init__() + + # Setup Inputs + self.inputs.append('x', '', 'X Variable') + self.inputs.append('y', '', 'Y Variable') + + # Setup Outputs + self.outputs.append('u', '', 'U Variable') + self.outputs.append('v', '', 'V Variable') + + # Simple model description + self.description = 'This model is a pass through)' + + self.availableDerivative = 1 + + # standard function call is y(, dydx, ...) = self.BlackBox(**{'x1':x1, 'x2':x2, ...}) + def BlackBox(self, *args, **kwargs): + # Note: This function as written only operates on the first input on a multi-input list + runCases, extras = self.parseInputs(*args, **kwargs) + x = [ + self.pyomo_value(runCases[i]['x']) for i in range(0, len(runCases)) + ] + y = [ + self.pyomo_value(runCases[i]['y']) for i in range(0, len(runCases)) + ] + + u = [] + dudx = [] + dudy = [] + v = [] + dvdx = [] + dvdy = [] + + for xval in x: + u.append(xval * units.dimensionless) + dudx.append(1.0 * units.dimensionless) + dudy.append(0.0 * units.dimensionless) + + for yval in y: + v.append(yval * units.dimensionless) + dvdx.append(0.0 * units.dimensionless) + dvdy.append(1.0 * units.dimensionless) + + if len(runCases) == 1: + returnMode = -2 + else: + returnMode = 1 + + if returnMode < 0: + returnMode = -1 * (returnMode + 1) + if returnMode == 0: + return [u[0], v[0]] + if returnMode == 1: + return [u[0], v[0]], [[dudx[0], dudy[0]], [dvdx[0], dvdy[0]]] + else: + if returnMode == 0: + opt = [] + for i in range(0, len(y)): + opt.append([u[i], v[i]]) + return opt + if returnMode == 1: + opt = [] + for i in range(0, len(y)): + opt.append( + [[u[i], v[i]], [[dudx[i], dudy[i]], [dvdx[i], dvdy[i]]]] + ) + return opt + + bb = PassThrough() + bbo = bb.BlackBox(1.0, 1.5) + bbs = bb.BlackBox_Standardized(1.0, 1.5) + + self.assertAlmostEqual(pyo.value(bbo[0][0]), pyo.value(bbs[0][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0][0]), pyomo_units.get_units(bbs[0][0]) + ) + self.assertAlmostEqual(pyo.value(bbo[0][1]), pyo.value(bbs[0][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0][1]), pyomo_units.get_units(bbs[0][1]) + ) + + self.assertAlmostEqual(pyo.value(bbo[1][0][0]), pyo.value(bbs[1][0][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][0]), pyomo_units.get_units(bbs[1][0][0]) + ) + self.assertAlmostEqual(pyo.value(bbo[1][0][1]), pyo.value(bbs[1][0][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][1]), pyomo_units.get_units(bbs[1][0][1]) + ) + self.assertAlmostEqual(pyo.value(bbo[1][1][0]), pyo.value(bbs[1][1][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][1][0]), pyomo_units.get_units(bbs[1][1][0]) + ) + self.assertAlmostEqual(pyo.value(bbo[1][1][1]), pyo.value(bbs[1][1][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][1][1]), pyomo_units.get_units(bbs[1][1][1]) + ) + + def test_edi_blackbox_standardized_2(self): + "Tests the blackbox_standardized function" + import pyomo + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel, ediprint + from pyomo.common.formatting import tostr + from pyomo.environ import units as pyomo_units + + class Norm_2(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the two norm' + self.inputs.append( + name='x', units='', description='The x variable', size=3 + ) + self.outputs.append(name='y', units='', description='The y variable') + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, *args, **kwargs): # The actual function that does things + runCases, remainingKwargs = self.parseInputs(*args, **kwargs) + + x = self.sanitizeInputs(runCases[0]['x']) + x = np.array([pyo.value(xval) for xval in x], dtype=np.float64) + + y = x[0] ** 2 + x[1] ** 2 + x[2] ** 2 # Compute y + dydx0 = 2 * x[0] # Compute dy/dx0 + dydx1 = 2 * x[1] # Compute dy/dx1 + dydx2 = 2 * x[2] # Compute dy/dx2 + + y = y * units.dimensionless + # dydx = np.array([dydx0,dydx1,dydx2]) * units.dimensionless + dydx0 = dydx0 * units.dimensionless + dydx1 = dydx1 * units.dimensionless + dydx2 = dydx2 * units.dimensionless + + return y, [[dydx0, dydx1, dydx2]] # return z, grad(z), hess(z)... + + bb = Norm_2() + bbo = bb.BlackBox(np.array([1, 1, 1])) + bbs = bb.BlackBox_Standardized(np.array([1, 1, 1])) + + self.assertAlmostEqual(pyo.value(bbs.values.y), pyo.value(bbo[0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbs.values.y), pyomo_units.get_units(bbo[0]) + ) + + self.assertAlmostEqual(pyo.value(bbs.first.y.x[0]), pyo.value(bbo[1][0][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbs.first.y.x[0]), pyomo_units.get_units(bbo[1][0][0]) + ) + + self.assertAlmostEqual(pyo.value(bbo[0]), pyo.value(bbs[0][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0]), pyomo_units.get_units(bbs[0][0]) + ) + + self.assertAlmostEqual(pyo.value(bbo[1][0][0]), pyo.value(bbs[1][0][0][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][0]), pyomo_units.get_units(bbs[1][0][0][0]) + ) + + self.assertAlmostEqual(pyo.value(bbo[1][0][1]), pyo.value(bbs[1][0][0][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][1]), pyomo_units.get_units(bbs[1][0][0][1]) + ) + + self.assertAlmostEqual(pyo.value(bbo[1][0][2]), pyo.value(bbs[1][0][0][2])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][2]), pyomo_units.get_units(bbs[1][0][0][2]) + ) + + def test_edi_blackbox_standardized_3(self): + "Tests the blackbox_standardized function" + import pyomo + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel, ediprint + from pyomo.common.formatting import tostr + from pyomo.environ import units as pyomo_units + + class VectorCast(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the two norm' + self.inputs.append(name='x', units='', description='The x variable') + self.outputs.append( + name='y', units='', description='The y variable', size=3 + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, *args, **kwargs): # The actual function that does things + runCases, remainingKwargs = self.parseInputs(*args, **kwargs) + + x = pyo.value(self.sanitizeInputs(runCases[0]['x'])) + + y = np.array([x, x, x]) * units.dimensionless + dydx = np.array([1.0, 1.0, 1.0]) * units.dimensionless + + return y, [dydx] # return z, grad(z), hess(z)... + + bb = VectorCast() + bbo = bb.BlackBox(2) + bbs = bb.BlackBox_Standardized(2) + + self.assertAlmostEqual(pyo.value(bbo[0][0]), pyo.value(bbs[0][0][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0][0]), pyomo_units.get_units(bbs[0][0][0]) + ) + + self.assertAlmostEqual(pyo.value(bbo[0][1]), pyo.value(bbs[0][0][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0][1]), pyomo_units.get_units(bbs[0][0][1]) + ) + + self.assertAlmostEqual(pyo.value(bbo[0][2]), pyo.value(bbs[0][0][2])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0][2]), pyomo_units.get_units(bbs[0][0][2]) + ) + + self.assertAlmostEqual(pyo.value(bbo[1][0][0]), pyo.value(bbs[1][0][0][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][0]), pyomo_units.get_units(bbs[1][0][0][0]) + ) + + self.assertAlmostEqual(pyo.value(bbo[1][0][1]), pyo.value(bbs[1][0][0][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][1]), pyomo_units.get_units(bbs[1][0][0][1]) + ) + + self.assertAlmostEqual(pyo.value(bbo[1][0][2]), pyo.value(bbs[1][0][0][2])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][2]), pyomo_units.get_units(bbs[1][0][0][2]) + ) + + def test_edi_blackbox_standardized_4(self): + "Tests the blackbox_standardized function" + import pyomo + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel, ediprint + from pyomo.common.formatting import tostr + from pyomo.environ import units as pyomo_units + + class Parabola(BlackBoxFunctionModel): + def __init__(self): + super(Parabola, self).__init__() + self.description = 'This model evaluates the function: y = x**2' + self.inputs.append( + name='x', size=3, units='', description='The x variable' + ) + self.outputs.append( + name='y', size=3, units='', description='The y variable' + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, *args, **kwargs): # The actual function that does things + runCases, remainingKwargs = self.parseInputs(*args, **kwargs) + + x = self.sanitizeInputs(runCases[0]['x']) + x = np.array([pyo.value(xval) for xval in x], dtype=np.float64) + + y = x**2 # Compute y + dydx = 2 * x # Compute dy/dx + + y = y + dydx = np.diag(dydx) + dydx = dydx + + return y, [dydx] # return z, grad(z), hess(z)... + + bb = Parabola() + bbo = bb.BlackBox(np.array([1, 1, 1])) + bbs = bb.BlackBox_Standardized(np.array([1, 1, 1])) + + self.assertAlmostEqual(pyo.value(bbo[0][0]), pyo.value(bbs[0][0][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0][0]), pyomo_units.get_units(bbs[0][0][0]) + ) + + self.assertAlmostEqual(pyo.value(bbo[0][1]), pyo.value(bbs[0][0][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0][1]), pyomo_units.get_units(bbs[0][0][1]) + ) + + self.assertAlmostEqual(pyo.value(bbo[0][2]), pyo.value(bbs[0][0][2])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0][2]), pyomo_units.get_units(bbs[0][0][2]) + ) + + self.assertAlmostEqual( + pyo.value(bbo[1][0][0][0]), pyo.value(bbs[1][0][0][0][0]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][0][0]), + pyomo_units.get_units(bbs[1][0][0][0][0]), + ) + + self.assertAlmostEqual( + pyo.value(bbo[1][0][0][1]), pyo.value(bbs[1][0][0][0][1]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][0][1]), + pyomo_units.get_units(bbs[1][0][0][0][1]), + ) + + self.assertAlmostEqual( + pyo.value(bbo[1][0][0][2]), pyo.value(bbs[1][0][0][0][2]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][0][2]), + pyomo_units.get_units(bbs[1][0][0][0][2]), + ) + + self.assertAlmostEqual( + pyo.value(bbo[1][0][1][0]), pyo.value(bbs[1][0][0][1][0]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][1][0]), + pyomo_units.get_units(bbs[1][0][0][1][0]), + ) + + self.assertAlmostEqual( + pyo.value(bbo[1][0][1][1]), pyo.value(bbs[1][0][0][1][1]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][1][1]), + pyomo_units.get_units(bbs[1][0][0][1][1]), + ) + + self.assertAlmostEqual( + pyo.value(bbo[1][0][1][2]), pyo.value(bbs[1][0][0][1][2]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][1][2]), + pyomo_units.get_units(bbs[1][0][0][1][2]), + ) + + self.assertAlmostEqual( + pyo.value(bbo[1][0][2][0]), pyo.value(bbs[1][0][0][2][0]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][2][0]), + pyomo_units.get_units(bbs[1][0][0][2][0]), + ) + + self.assertAlmostEqual( + pyo.value(bbo[1][0][2][1]), pyo.value(bbs[1][0][0][2][1]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][2][1]), + pyomo_units.get_units(bbs[1][0][0][2][1]), + ) + + self.assertAlmostEqual( + pyo.value(bbo[1][0][2][2]), pyo.value(bbs[1][0][0][2][2]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][2][2]), + pyomo_units.get_units(bbs[1][0][0][2][2]), + ) + + def test_edi_blackbox_standardized_5(self): + "Tests the blackbox_standardized function" + import pyomo + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel, ediprint + from pyomo.common.formatting import tostr + from pyomo.environ import units as pyomo_units + + class MatrixSum(BlackBoxFunctionModel): + def __init__(self): + # Set up all the attributes by calling Model.__init__ + super().__init__() + + # Setup Inputs + self.inputs.append('x', '', 'Independent Variable', size=[2, 2]) + + # Setup Outputs + self.outputs.append('y', '', 'Dependent Variable') + + # Simple model description + self.description = 'This model sums all elements of a 2x2 matrix' + + self.availableDerivative = 1 + + # standard function call is y(, dydx, ...) = self.BlackBox(**{'x1':x1, 'x2':x2, ...}) + def BlackBox(self, *args, **kwargs): + runCases, extras = self.parseInputs(*args, **kwargs) + x = [ + self.pyomo_value(runCases[i]['x']) for i in range(0, len(runCases)) + ] + + y = np.sum(x) * units.dimensionless + dydx = np.ones([2, 2]) * units.dimensionless + + return y, dydx + + bb = MatrixSum() + bbo = bb.BlackBox(np.array([[1, 1], [1, 1]])) + bbs = bb.BlackBox_Standardized(np.array([[1, 1], [1, 1]])) + + self.assertAlmostEqual(pyo.value(bbo[0]), pyo.value(bbs[0][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0]), pyomo_units.get_units(bbs[0][0]) + ) + + self.assertAlmostEqual(pyo.value(bbo[1][0][0]), pyo.value(bbs[1][0][0][0][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][0]), + pyomo_units.get_units(bbs[1][0][0][0][0]), + ) + + self.assertAlmostEqual(pyo.value(bbo[1][0][1]), pyo.value(bbs[1][0][0][0][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][1]), + pyomo_units.get_units(bbs[1][0][0][0][1]), + ) + + self.assertAlmostEqual(pyo.value(bbo[1][1][0]), pyo.value(bbs[1][0][0][1][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][1][0]), + pyomo_units.get_units(bbs[1][0][0][1][0]), + ) + + self.assertAlmostEqual(pyo.value(bbo[1][1][1]), pyo.value(bbs[1][0][0][1][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][1][1]), + pyomo_units.get_units(bbs[1][0][0][1][1]), + ) + + def test_edi_blackbox_standardized_6(self): + "Tests the blackbox_standardized function" + import pyomo + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel, ediprint + from pyomo.common.formatting import tostr + from pyomo.environ import units as pyomo_units + + class MatrixCast(BlackBoxFunctionModel): + def __init__(self): + # Set up all the attributes by calling Model.__init__ + super().__init__() + + # Setup Inputs + self.inputs.append('x', '', 'Independent Variable') + + # Setup Outputs + self.outputs.append('y', '', 'Dependent Variable', size=[2, 2]) + + # Simple model description + self.description = 'This model casts a number to a 2x2 matrix' + + self.availableDerivative = 1 + + # standard function call is y(, dydx, ...) = self.BlackBox(**{'x1':x1, 'x2':x2, ...}) + def BlackBox(self, *args, **kwargs): + runCases, extras = self.parseInputs(*args, **kwargs) + x = [ + self.pyomo_value(runCases[i]['x']) for i in range(0, len(runCases)) + ] + + y = x[0] * np.ones([2, 2]) * units.dimensionless + dydx = np.ones([2, 2]) * units.dimensionless + + return y, dydx + + bb = MatrixCast() + bbo = bb.BlackBox(2) + bbs = bb.BlackBox_Standardized(2) + + self.assertAlmostEqual(pyo.value(bbo[0][0][0]), pyo.value(bbs[0][0][0][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0][0][0]), pyomo_units.get_units(bbs[0][0][0][0]) + ) + + self.assertAlmostEqual(pyo.value(bbo[0][0][1]), pyo.value(bbs[0][0][0][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0][0][1]), pyomo_units.get_units(bbs[0][0][0][1]) + ) + + self.assertAlmostEqual(pyo.value(bbo[0][1][0]), pyo.value(bbs[0][0][1][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0][1][0]), pyomo_units.get_units(bbs[0][0][1][0]) + ) + + self.assertAlmostEqual(pyo.value(bbo[0][1][1]), pyo.value(bbs[0][0][1][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[0][1][1]), pyomo_units.get_units(bbs[0][0][1][1]) + ) + + self.assertAlmostEqual(pyo.value(bbo[1][0][0]), pyo.value(bbs[1][0][0][0][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][0]), + pyomo_units.get_units(bbs[1][0][0][0][0]), + ) + + self.assertAlmostEqual(pyo.value(bbo[1][0][1]), pyo.value(bbs[1][0][0][0][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][0][1]), + pyomo_units.get_units(bbs[1][0][0][0][1]), + ) + + self.assertAlmostEqual(pyo.value(bbo[1][1][0]), pyo.value(bbs[1][0][0][1][0])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][1][0]), + pyomo_units.get_units(bbs[1][0][0][1][0]), + ) + + self.assertAlmostEqual(pyo.value(bbo[1][1][1]), pyo.value(bbs[1][0][0][1][1])) + self.assertAlmostEqual( + pyomo_units.get_units(bbo[1][1][1]), + pyomo_units.get_units(bbs[1][0][0][1][1]), + ) + + def test_edi_blackbox_multicase(self): + "Tests the multicase function" + import pyomo + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel, ediprint + from pyomo.common.formatting import tostr + from pyomo.environ import units as pyomo_units + + class VectorCast(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the two norm' + self.inputs.append(name='x', units='', description='The x variable') + self.outputs.append( + name='y', units='', description='The y variable', size=3 + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, *args, **kwargs): # The actual function that does things + runCases, remainingKwargs = self.parseInputs(*args, **kwargs) + + x = pyo.value(self.sanitizeInputs(runCases[0]['x'])) + + y = np.array([x, x, x]) * units.dimensionless + dydx = np.array([1.0, 1.0, 1.0]) * units.dimensionless + + return y, [dydx] # return z, grad(z), hess(z)... + + bb = VectorCast() + bbs = bb.BlackBox_Standardized([[1], [2]]) + bbm = bb.MultiCase([[1], [2]]) + + self.assertAlmostEqual( + pyo.value(bbs[0].values.y[0]), pyo.value(bbm[0].values.y[0]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbs[0].values.y[0]), + pyomo_units.get_units(bbm[0].values.y[0]), + ) + + self.assertAlmostEqual( + pyo.value(bbs[0].values.y[1]), pyo.value(bbm[0].values.y[1]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbs[0].values.y[1]), + pyomo_units.get_units(bbm[0].values.y[1]), + ) + + self.assertAlmostEqual( + pyo.value(bbs[0].values.y[2]), pyo.value(bbm[0].values.y[2]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbs[0].values.y[2]), + pyomo_units.get_units(bbm[0].values.y[2]), + ) + + self.assertAlmostEqual( + pyo.value(bbs[0].first.y.x[0]), pyo.value(bbm[0].first.y.x[0]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbs[0].first.y.x[0]), + pyomo_units.get_units(bbm[0].first.y.x[0]), + ) + + self.assertAlmostEqual( + pyo.value(bbs[0].first.y.x[1]), pyo.value(bbm[0].first.y.x[1]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbs[0].first.y.x[1]), + pyomo_units.get_units(bbm[0].first.y.x[1]), + ) + + self.assertAlmostEqual( + pyo.value(bbs[0].first.y.x[2]), pyo.value(bbm[0].first.y.x[2]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbs[0].first.y.x[2]), + pyomo_units.get_units(bbm[0].first.y.x[2]), + ) + + self.assertAlmostEqual( + pyo.value(bbs[1].values.y[0]), pyo.value(bbm[1].values.y[0]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbs[1].values.y[0]), + pyomo_units.get_units(bbm[1].values.y[0]), + ) + + self.assertAlmostEqual( + pyo.value(bbs[1].values.y[1]), pyo.value(bbm[1].values.y[1]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbs[1].values.y[1]), + pyomo_units.get_units(bbm[1].values.y[1]), + ) + + self.assertAlmostEqual( + pyo.value(bbs[1].values.y[2]), pyo.value(bbm[1].values.y[2]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbs[1].values.y[2]), + pyomo_units.get_units(bbm[1].values.y[2]), + ) + + self.assertAlmostEqual( + pyo.value(bbs[1].first.y.x[0]), pyo.value(bbm[1].first.y.x[0]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbs[1].first.y.x[0]), + pyomo_units.get_units(bbm[1].first.y.x[0]), + ) + + self.assertAlmostEqual( + pyo.value(bbs[1].first.y.x[1]), pyo.value(bbm[1].first.y.x[1]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbs[1].first.y.x[1]), + pyomo_units.get_units(bbm[1].first.y.x[1]), + ) + + self.assertAlmostEqual( + pyo.value(bbs[1].first.y.x[2]), pyo.value(bbm[1].first.y.x[2]) + ) + self.assertAlmostEqual( + pyomo_units.get_units(bbs[1].first.y.x[2]), + pyomo_units.get_units(bbm[1].first.y.x[2]), + ) + + +if __name__ == '__main__': + unittest.main() diff --git a/pyomo/contrib/edi/tests/test_docSnippets.py b/pyomo/contrib/edi/tests/test_docSnippets.py new file mode 100644 index 00000000000..09bc1fa3ab4 --- /dev/null +++ b/pyomo/contrib/edi/tests/test_docSnippets.py @@ -0,0 +1,1161 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# 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. +# +# Development of this module was conducted as part of the Institute for +# the Design of Advanced Energy Systems (IDAES) with support through the +# Simulation-Based Engineering, Crosscutting Research Program within the +# U.S. Department of Energy’s Office of Fossil Energy and Carbon Management. +# +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import pyomo.common.unittest as unittest +import pyomo.environ as pyo +from pyomo.common.dependencies import attempt_import + +from pyomo.core.base.units_container import pint_available + +from pyomo.common.dependencies import numpy, numpy_available +from pyomo.common.dependencies import scipy, scipy_available + +egb, egb_available = attempt_import( + "pyomo.contrib.pynumero.interfaces.external_grey_box" +) + +formulation_available = False +try: + from pyomo.contrib.edi import Formulation + + formulation_available = True +except: + pass + # formulation_available = False + +blackbox_available = False +try: + from pyomo.contrib.edi import BlackBoxFunctionModel + + blackbox_available = True +except: + pass + # blackbox_available = False + +if numpy_available: + import numpy as np + + +@unittest.skipIf( + not egb_available, 'Testing pyomo.contrib.edi requires pynumero external grey boxes' +) +@unittest.skipIf(not formulation_available, 'Formulation import failed') +@unittest.skipIf(not blackbox_available, 'Blackbox import failed') +@unittest.skipIf(not numpy_available, 'Testing pyomo.contrib.edi requires numpy') +@unittest.skipIf(not scipy_available, 'Testing pyomo.contrib.edi requires scipy') +@unittest.skipIf(not pint_available, 'Testing units requires pint') +class TestEDISnippets(unittest.TestCase): + def test_edi_snippet_formuation_01(self): + # BEGIN: Formulation_Snippet_01 + from pyomo.contrib.edi import Formulation + + f = Formulation() + # END: Formulation_Snippet_01 + + def test_edi_snippet_formuation_02(self): + # BEGIN: Formulation_Snippet_02 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='') + # END: Formulation_Snippet_02 + + def test_edi_snippet_formuation_03(self): + # BEGIN: Formulation_Snippet_03 + from pyomo.contrib.edi import Formulation + + f = Formulation() + c = f.Constant(name='c', value=1.0, units='') + # END: Formulation_Snippet_03 + + def test_edi_snippet_formuation_04(self): + # BEGIN: Formulation_Snippet_04 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='') + y = f.Variable(name='y', guess=1.0, units='') + c = f.Constant(name='c', value=1.0, units='') + f.Objective(c * x + y) + # END: Formulation_Snippet_04 + + def test_edi_snippet_formuation_05(self): + # BEGIN: Formulation_Snippet_05 + from pyomo.contrib.edi import Formulation + from pyomo.environ import maximize, minimize + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='') + y = f.Variable(name='y', guess=1.0, units='') + c = f.Constant(name='c', value=1.0, units='') + f.Objective(c * x + y, sense=maximize) + # END: Formulation_Snippet_05 + + def test_edi_snippet_formuation_06(self): + # BEGIN: Formulation_Snippet_06 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='') + y = f.Variable(name='y', guess=1.0, units='') + c = f.Constant(name='c', value=1.0, units='') + f.Objective(c * x + y) + f.Constraint(x**2 + y**2 <= 1.0) + f.Constraint(x >= 0) + f.Constraint(y <= 0) + # END: Formulation_Snippet_06 + + def test_edi_snippet_formuation_07(self): + # BEGIN: Formulation_Snippet_07 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='') + y = f.Variable(name='y', guess=1.0, units='') + c = f.Constant(name='c', value=1.0, units='') + f.Objective(c * x + y) + f.ConstraintList([x**2 + y**2 <= 1.0, x >= 0, y <= 0]) + # END: Formulation_Snippet_07 + + def test_edi_snippet_formuation_08(self): + # BEGIN: Formulation_Snippet_08 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='') + y = f.Variable(name='y', guess=1.0, units='') + c = f.Constant(name='c', value=1.0, units='') + f.Objective(c * x + y) + + constraintList = [x**2 + y**2 <= 1.0, x >= 0, y <= 0] + + f.ConstraintList(constraintList) + # END: Formulation_Snippet_08 + + def test_edi_snippet_formuation_09(self): + # BEGIN: Formulation_Snippet_09 + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='y variable') + z = f.Variable(name='z', guess=1.0, units='m^2', description='Output var') + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', units='ft**2', description='Output variable' + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + # Converts to correct units then casts to float + x = pyo.value(units.convert(x, self.inputs['x'].units)) + y = pyo.value(units.convert(y, self.inputs['y'].units)) + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + f.Constraint(z <= 1 * units.m**2) + + f.RuntimeConstraint(z, '==', [x, y], UnitCircle()) + # END: Formulation_Snippet_09 + + def test_edi_snippet_formuation_14(self): + # BEGIN: Formulation_Snippet_14 + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='y variable') + z = f.Variable(name='z', guess=1.0, units='m^2', description='Output var') + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', units='ft**2', description='Output variable' + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + # Converts to correct units then casts to float + x = pyo.value(units.convert(x, self.inputs[0].units)) + y = pyo.value(units.convert(y, self.inputs[1].units)) + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + f.ConstraintList([z <= 1 * units.m**2, [z, '==', [x, y], UnitCircle()]]) + # END: Formulation_Snippet_14 + + ### This will fail validatation, but should construct appropriately + + # BEGIN: Formulation_Snippet_15 + f.ConstraintList([z <= 1 * units.m**2, (z, '==', [x, y], UnitCircle())]) + # END: Formulation_Snippet_15 + + # BEGIN: Formulation_Snippet_16 + f.ConstraintList([z <= 1 * units.m**2, [z, '==', [x, y], UnitCircle()]]) + # END: Formulation_Snippet_16 + + # BEGIN: Formulation_Snippet_17 + f.ConstraintList( + [ + z <= 1 * units.m**2, + { + 'outputs': z, + 'operators': '==', + 'inputs': [x, y], + 'black_box': UnitCircle(), + }, + ] + ) + # END: Formulation_Snippet_17 + + # BEGIN: Formulation_Snippet_18 + f.ConstraintList([z <= 1 * units.m**2, ([z], ['=='], [x, y], UnitCircle())]) + # END: Formulation_Snippet_18 + + def test_edi_snippet_variables_01(self): + # BEGIN: Variables_Snippet_01 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + # END: Variables_Snippet_01 + + def test_edi_snippet_variables_02(self): + # BEGIN: Variables_Snippet_02 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable('x', 1.0, 'm') + # END: Variables_Snippet_02 + + def test_edi_snippet_variables_03(self): + # BEGIN: Variables_Snippet_03 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable( + name='x', + guess=1.0, + units='m', + description='The x variable', + bounds=[-10, 10], + ) + # END: Variables_Snippet_03 + + def test_edi_snippet_variables_04(self): + # BEGIN: Variables_Snippet_04 + from pyomo.contrib.edi import Formulation + from pyomo.environ import Integers + + f = Formulation() + x = f.Variable( + name='x', + guess=1.0, + units='m', + description='The x variable', + domain=Integers, + ) + # END: Variables_Snippet_04 + + def test_edi_snippet_variables_05(self): + # BEGIN: Variables_Snippet_05 + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units=units.m, description='The x variable') + # END: Variables_Snippet_05 + + def test_edi_snippet_variables_06(self): + # BEGIN: Variables_Snippet_06 + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable( + name='x', guess=1.0, units='m', description='The x variable', size=5 + ) + # END: Variables_Snippet_06 + + def test_edi_snippet_variables_07(self): + # BEGIN: Variables_Snippet_07 + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable( + name='x', guess=1.0, units='m', description='The x variable', size=[10, 2] + ) + # END: Variables_Snippet_07 + + def test_edi_snippet_variables_08(self): + # BEGIN: Variables_Snippet_08 + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable( + name='x', guess=1.0, units='kg*m/s**2', description='The x variable' + ) + # END: Variables_Snippet_08 + + def test_edi_snippet_constants_01(self): + # BEGIN: Constants_Snippet_01 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Constant(name='c', value=1.0, units='m', description='A constant c') + # END: Constants_Snippet_01 + + def test_edi_snippet_constants_02(self): + # BEGIN: Constants_Snippet_02 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Constant('c', 1.0, 'm') + # END: Constants_Snippet_02 + + def test_edi_snippet_constants_03(self): + # BEGIN: Constants_Snippet_03 + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Constant(name='c', value=1.0, units=units.m, description='A constant c') + # END: Constants_Snippet_03 + + def test_edi_snippet_constants_04(self): + # BEGIN: Constants_Snippet_04 + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Constant( + name='c', value=1.0, units='m', description='A constant c', size=5 + ) + # END: Constants_Snippet_04 + + def test_edi_snippet_constants_05(self): + # BEGIN: Constants_Snippet_05 + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Constant( + name='c', value=1.0, units='m', description='A constant c', size=[10, 2] + ) + # END: Constants_Snippet_05 + + def test_edi_snippet_constants_06(self): + # BEGIN: Constants_Snippet_06 + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Constant( + name='c', value=1.0, units='kg*m/s**2', description='A constant c' + ) + # END: Constants_Snippet_06 + + def test_edi_snippet_objectives_01(self): + # BEGIN: Objectives_Snippet_01 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + c = f.Constant(name='c', value=1.0, units='', description='A constant c') + f.Objective(c * x + y) # Default is minimize + # END: Objectives_Snippet_01 + + def test_edi_snippet_objectives_02(self): + # BEGIN: Objectives_Snippet_02 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + c = f.Constant(name='c', value=1.0, units='', description='A constant c') + f.Objective(c * x**4 + y**4) # Default is minimize + # END: Objectives_Snippet_02 + + def test_edi_snippet_objectives_03(self): + # BEGIN: Objectives_Snippet_03 + from pyomo.contrib.edi import Formulation + from pyomo.environ import minimize, maximize + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + c = f.Constant(name='c', value=1.0, units='', description='A constant c') + f.Objective(c * x**4 + y**4, sense=minimize) + # END: Objectives_Snippet_03 + + def test_edi_snippet_objectives_04(self): + # BEGIN: Objectives_Snippet_04 + from pyomo.contrib.edi import Formulation + from pyomo.environ import minimize, maximize + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + c = f.Constant(name='c', value=1.0, units='', description='A constant c') + f.Objective(c * x**4 + y**4, sense=1) # 1 corresponds to minimize + # END: Objectives_Snippet_04 + + def test_edi_snippet_objectives_05(self): + # BEGIN: Objectives_Snippet_05 + from pyomo.contrib.edi import Formulation + from pyomo.environ import minimize, maximize + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + c = f.Constant(name='c', value=1.0, units='', description='A constant c') + f.Objective(-c * x**4 - y**4, sense=maximize) + # END: Objectives_Snippet_05 + + def test_edi_snippet_objectives_06(self): + # BEGIN: Objectives_Snippet_06 + from pyomo.contrib.edi import Formulation + from pyomo.environ import minimize, maximize + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + c = f.Constant(name='c', value=1.0, units='', description='A constant c') + f.Objective(-c * x**4 - y**4, sense=-1) # -1 corresponds to maximize + # END: Objectives_Snippet_06 + + def test_edi_snippet_objectives_07(self): + # BEGIN: Objectives_Snippet_07 + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable( + name='x', + guess=1.0, + units='m', + description='The x variable', + bounds=[0, 100], + size=3, + ) + y = f.Variable( + name='y', guess=1.0, units='m', description='The y variable', size=[2, 2] + ) + c = f.Constant( + name='c', value=1.0, units='', description='A constant c', size=3 + ) + f.Objective( + c[0] * x[0] + + c[1] * x[1] + + c[2] * x[2] + + y[0, 0] ** 4 + + y[0, 1] ** 4 + + y[1, 0] ** 4 + + y[1, 1] ** 4 + ) # Default is minimize + # END: Objectives_Snippet_07 + + def test_edi_snippet_constraints_01(self): + # BEGIN: Constraints_Snippet_01 + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + c = f.Constant(name='c', value=1.0, units='', description='A constant c') + f.Objective(c * x + y) + f.ConstraintList([x**2 + y**2 <= 1.0 * units.m**2, x <= 0.75 * units.m, x >= y]) + # END: Constraints_Snippet_01 + + def test_edi_snippet_constraints_02(self): + # BEGIN: Constraints_Snippet_02 + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + c = f.Constant(name='c', value=1.0, units='', description='A constant c') + f.Objective(c * x + y) + f.Constraint(x**2 + y**2 <= 1.0 * units.m**2) + f.Constraint(x <= 0.75 * units.m) + f.Constraint(x >= y) + # END: Constraints_Snippet_02 + + def test_edi_snippet_constraints_03(self): + # BEGIN: Constraints_Snippet_03 + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable( + name='x', + guess=1.0, + units='m', + description='The x variable', + bounds=[0, 100], + size=3, + ) + y = f.Variable( + name='y', guess=1.0, units='m', description='The y variable', size=[2, 2] + ) + c = f.Constant( + name='c', value=1.0, units='', description='A constant c', size=3 + ) + f.Objective( + c[0] * x[0] + + c[1] * x[1] + + c[2] * x[2] + + y[0, 0] ** 4 + + y[0, 1] ** 4 + + y[1, 0] ** 4 + + y[1, 1] ** 4 + ) # Default is minimize + f.ConstraintList( + [ + x[0] ** 2 + x[1] ** 2 + x[2] ** 2 <= 1.0 * units.m, + y[0, 0] >= 1.0 * units.m, + y[0, 1] >= 1.0 * units.m, + y[1, 0] >= 1.0 * units.m, + y[1, 1] >= 1.0 * units.m, + x[0] >= y[0, 0], + ] + ) + # END: Constraints_Snippet_03 + + def test_edi_snippet_runtimeconstraints_01(self): + # BEGIN: RuntimeConstraints_Snippet_01 + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import BlackBoxFunctionModel + + class Parabola(BlackBoxFunctionModel): + def __init__(self): + # Call parent init + super().__init__() + + # A brief description of the model + self.description = 'This model evaluates the function: y = x**2' + + # Append the model inputs + self.inputs.append(name='x', units='ft', description='The x variable') + + # Append the model outputs + self.outputs.append( + name='y', units='ft**2', description='The y variable' + ) + + # Set the highest available derivative + # Should be 1 for most cases but defaults to 0 + self.availableDerivative = 1 + + def BlackBox(self, x): # The actual function that does things + # Convert to correct units and cast to a float + x = pyo.value(units.convert(x, self.inputs['x'].units)) + + # Compute y + y = x**2 + + # Compute dy/dx + dydx = 2 * x + + # Add the units to the output + y = y * self.outputs['y'].units + + # Add the units to the derivative for output + dydx = dydx * self.outputs['y'].units / self.inputs['x'].units + + # Return using the output packing guidelines described in the documentation: + # returnVal[0] = output_value + # returnVal[1] = jacobian + # returnVal[1][0] = derivative_scalarOutput_wrt_0th_input + return y, [dydx] + + # END: RuntimeConstraints_Snippet_01 + + def test_edi_snippet_runtimeconstraints_02(self): + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import BlackBoxFunctionModel + + class Parabola(BlackBoxFunctionModel): + def __init__(self): + # BEGIN: RuntimeConstraints_Snippet_02 + # Call parent init + super().__init__() + # END: RuntimeConstraints_Snippet_02 + + # A brief description of the model + self.description = 'This model evaluates the function: y = x**2' + + # Append the model inputs + self.inputs.append(name='x', units='ft', description='The x variable') + + # Append the model outputs + self.outputs.append( + name='y', units='ft**2', description='The y variable' + ) + + # BEGIN: RuntimeConstraints_Snippet_05 + # Set the highest available derivative + # Should be 1 for most cases but defaults to 0 + self.availableDerivative = 1 + # END: RuntimeConstraints_Snippet_05 + + def BlackBox(self, x): # The actual function that does things + storeX = x + + # BEGIN: RuntimeConstraints_Snippet_06 + x = units.convert(x, self.inputs['x'].units) + # END: RuntimeConstraints_Snippet_06 + + x = storeX + + # BEGIN: RuntimeConstraints_Snippet_07 + # Convert to correct units and cast to a float + x = pyo.value(units.convert(x, self.inputs['x'].units)) + # END: RuntimeConstraints_Snippet_07 + + # Compute y + y = x**2 + + # Compute dy/dx + dydx = 2 * x + + # BEGIN: RuntimeConstraints_Snippet_08 + # Add the units to the output + y = y * self.outputs['y'].units + + # Add the units to the derivative for output + dydx = dydx * self.outputs['y'].units / self.inputs['x'].units + # END: RuntimeConstraints_Snippet_08 + + # BEGIN: RuntimeConstraints_Snippet_09 + # Return using the output packing guidelines described in the documentation: + # returnVal[0] = output_value + # returnVal[1] = jacobian + # returnVal[1][0] = derivative_scalarOutput_wrt_0th_input + return y, [dydx] + # END: RuntimeConstraints_Snippet_09 + + def test_edi_snippet_runtimeconstraints_03(self): + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + z = f.Variable(name='z', guess=1.0, units='m^2', description='The z variable') + + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + # BEGIN: RuntimeConstraints_Snippet_03 + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + # END: RuntimeConstraints_Snippet_03 + self.outputs.append( + name='z', + units='ft**2', + description='Resultant of the unit circle evaluation', + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + x = pyo.value(units.convert(x, self.inputs['x'].units)) + y = pyo.value(units.convert(y, self.inputs['y'].units)) + z = x**2 + y**2 + dzdx = 2 * x + dzdy = 2 * y + z = z * self.outputs['z'].units + dzdx = dzdx * self.outputs['z'].units / self.inputs['x'].units + dzdy = dzdy * self.outputs['z'].units / self.inputs['y'].units + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + f.ConstraintList([(z, '==', [x, y], UnitCircle()), z <= 1 * units.m**2]) + + def test_edi_snippet_runtimeconstraints_04(self): + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + from pyomo.common.formatting import tostr + + class PassThrough(BlackBoxFunctionModel): + def __init__(self): + # Set up all the attributes by calling Model.__init__ + super().__init__() + + # Setup Inputs + self.inputs.append('x', '', 'X Variable') + self.inputs.append('y', '', 'Y Variable') + + # BEGIN: RuntimeConstraints_Snippet_04 + # Setup Outputs + self.outputs.append('u', '', 'U Variable') + self.outputs.append('v', '', 'V Variable') + # END: RuntimeConstraints_Snippet_04 + + # Simple model description + self.description = 'This model is a pass through)' + + self.availableDerivative = 1 + + # standard function call is y(, dydx, ...) = self.BlackBox(**{'x1':x1, 'x2':x2, ...}) + def BlackBox(self, *args, **kwargs): + runCases, extras = self.parseInputs(*args, **kwargs) + x = [ + self.pyomo_value(runCases[i]['x']) for i in range(0, len(runCases)) + ] + y = [ + self.pyomo_value(runCases[i]['x']) for i in range(0, len(runCases)) + ] + + u = [] + dudx = [] + dudy = [] + v = [] + dvdx = [] + dvdy = [] + + for xval in x: + u.append(xval * units.dimensionless) + dudx.append(1.0 * units.dimensionless) + dudy.append(0.0 * units.dimensionless) + + for yval in y: + v.append(yval * units.dimensionless) + dvdx.append(0.0 * units.dimensionless) + dvdy.append(1.0 * units.dimensionless) + + returnMode = 1 + + if returnMode < 0: + returnMode = -1 * (returnMode + 1) + if returnMode == 0: + return [u[0], v[0]] + if returnMode == 1: + return [u[0], v[0]], [[dudx[0], dudy[0]], [dvdx[0], dvdy[0]]] + else: + if returnMode == 0: + opt = [] + for i in range(0, len(y)): + opt.append([u[i], v[i]]) + return opt + if returnMode == 1: + opt = [] + for i in range(0, len(y)): + opt.append( + [[u[i], v[i]], [[dudx[i], dudy[i]], [dvdx[i], dvdy[i]]]] + ) + return opt + + bb = PassThrough() + bbo = bb.BlackBox(1.0, 1.0) + bbo = bb.BlackBox({'x': np.linspace(0, 10, 11), 'y': np.linspace(0, 10, 11)}) + + def test_edi_snippet_runtimeconstraints_10(self): + # BEGIN: RuntimeConstraints_Snippet_10 + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + z = f.Variable(name='z', guess=1.0, units='m^2', description='The z variable') + + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', + units='ft**2', + description='Resultant of the unit circle evaluation', + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + x = pyo.value(units.convert(x, self.inputs['x'].units)) + y = pyo.value(units.convert(y, self.inputs['y'].units)) + z = x**2 + y**2 + dzdx = 2 * x + dzdy = 2 * y + z = z * self.outputs['z'].units + dzdx = dzdx * self.outputs['z'].units / self.inputs['x'].units + dzdy = dzdy * self.outputs['z'].units / self.inputs['y'].units + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + f.ConstraintList([(z, '==', [x, y], UnitCircle()), z <= 1 * units.m**2]) + # END: RuntimeConstraints_Snippet_10 + + def test_edi_snippet_runtimeconstraints_11(self): + # BEGIN: RuntimeConstraints_Snippet_11 + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + z = f.Variable(name='z', guess=1.0, units='m^2', description='The z variable') + + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', + units='ft**2', + description='Resultant of the unit circle evaluation', + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + x = pyo.value(units.convert(x, self.inputs['x'].units)) + y = pyo.value(units.convert(y, self.inputs['y'].units)) + z = x**2 + y**2 + dzdx = 2 * x + dzdy = 2 * y + z = z * self.outputs['z'].units + dzdx = dzdx * self.outputs['z'].units / self.inputs['x'].units + dzdy = dzdy * self.outputs['z'].units / self.inputs['y'].units + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + m = UnitCircle() + bb = m.BlackBox(0.5 * units.m, 0.5 * units.m) + bbs = m.BlackBox_Standardized(0.5 * units.m, 0.5 * units.m) + # END: RuntimeConstraints_Snippet_11 + + def test_edi_snippet_runtimeconstraints_12(self): + # BEGIN: RuntimeConstraints_Snippet_12 + from collections import namedtuple + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', + units='ft**2', + description='Resultant of the unit circle evaluation', + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + x = pyo.value(units.convert(x, self.inputs['x'].units)) + y = pyo.value(units.convert(y, self.inputs['y'].units)) + z = x**2 + y**2 + dzdx = 2 * x + dzdy = 2 * y + z = z * self.outputs['z'].units + dzdx = dzdx * self.outputs['z'].units / self.inputs['x'].units + dzdy = dzdy * self.outputs['z'].units / self.inputs['y'].units + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + def BlackBox_Standardized(self, x, y): # The standardized function + res = self.BlackBox(x, y) + + returnTuple = namedtuple('returnTuple', ['values', 'first', 'second']) + optTuple = namedtuple('optTuple', [o.name for o in self.outputs]) + iptTuple = namedtuple('iptTuple', [i.name for i in self.inputs]) + + values = optTuple(res[0]) + first = optTuple(iptTuple(res[1][0], res[1][1])) + second = None # Second derivatives not currently supported + + return returnTuple(values, first, second) + + m = UnitCircle() + bb = m.BlackBox(0.5 * units.m, 0.5 * units.m) + bbs = m.BlackBox_Standardized(0.5 * units.m, 0.5 * units.m) + # END: RuntimeConstraints_Snippet_12 + + def test_edi_snippet_runtimeconstraints_13(self): + # BEGIN: RuntimeConstraints_Snippet_13 + from collections import namedtuple + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super().__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', + units='ft**2', + description='Resultant of the unit circle evaluation', + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + x = pyo.value(units.convert(x, self.inputs['x'].units)) + y = pyo.value(units.convert(y, self.inputs['y'].units)) + z = x**2 + y**2 + dzdx = 2 * x + dzdy = 2 * y + z = z * self.outputs['z'].units + dzdx = dzdx * self.outputs['z'].units / self.inputs['x'].units + dzdy = dzdy * self.outputs['z'].units / self.inputs['y'].units + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + def MultiCase(self, list_of_cases): # The multi case function + returnTuple = namedtuple('returnTuple', ['values', 'first', 'second']) + optTuple = namedtuple('optTuple', [o.name for o in self.outputs]) + iptTuple = namedtuple('iptTuple', [i.name for i in self.inputs]) + + outputList = [] + + for i in range(0, len(list_of_cases)): # Could parallelize this loop + res = self.BlackBox_Standardized( + list_of_cases[i][0], list_of_cases[i][1] + ) + values = optTuple(res.values) + first = optTuple(iptTuple(res.first.z.x, res.first.z.y)) + second = None # Second derivatives not currently supported + outputList.append(returnTuple(values, first, second)) + + return outputList + + m = UnitCircle() + bb = m.BlackBox(0.5 * units.m, 0.5 * units.m) + bbm = m.MultiCase( + [[0.5 * units.m, 0.5 * units.m], [1.0 * units.m, 1.0 * units.m]] + ) + # END: RuntimeConstraints_Snippet_13 + + def test_edi_snippet_advancedRTC_01(self): + # BEGIN: AdvancedRTC_Snippet_01 + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + class SignomialTest(BlackBoxFunctionModel): + def __init__(self): + # Set up all the attributes by calling Model.__init__ + super().__init__() + + # Setup Inputs + self.inputs.append('x', '', 'Independent Variable') + + # Setup Outputs + self.outputs.append('y', '', 'Dependent Variable') + + # Simple model description + self.description = ( + 'This model evaluates the ' + 'function: max([-6*x-6, x**4-3*x**2])' + ) + + self.availableDerivative = 1 + + def BlackBox(self, *args, **kwargs): + runCases, extras = self.parseInputs(*args, **kwargs) + + x = np.array( + [pyo.value(runCases[i]['x']) for i in range(0, len(runCases))] + ) + + y = np.maximum(-6 * x - 6, x**4 - 3 * x**2) + dydx = 4 * x**3 - 6 * x + ddy_ddx = 12 * x**2 - 6 + gradientSwitch = -6 * x - 6 > x**4 - 3 * x**2 + dydx[gradientSwitch] = -6 + ddy_ddx[gradientSwitch] = 0 + + y = [yval * units.dimensionless for yval in y] + dydx = [dydx[i] * units.dimensionless for i in range(0, len(dydx))] + + returnMode = 1 + + if returnMode < 0: + returnMode = -1 * (returnMode + 1) + if returnMode == 0: + return y[0] + if returnMode == 1: + return y[0], dydx + else: + if returnMode == 0: + opt = [] + for i in range(0, len(y)): + opt.append([y[i]]) + return opt + if returnMode == 1: + opt = [] + for i in range(0, len(y)): + opt.append([[y[i]], [[[dydx[i]]]]]) + return opt + + s = SignomialTest() + ivals = [[x] for x in np.linspace(-2, 2, 11)] + + # How the black box may be called using EDI + bbo = s.BlackBox(**{'x': 0.5}) + bbo = s.BlackBox({'x': 0.5}) + bbo = s.BlackBox(**{'x': 0.5, 'optn': True}) + + # Additional options available with parseInputs + bbo = s.BlackBox(*[0.5], **{'optn1': True, 'optn2': False}) + bbo = s.BlackBox(*[0.5, True], **{'optn': False}) + bbo = s.BlackBox({'x': np.linspace(-2, 2, 11)}) + bbo = s.BlackBox({'x': [x for x in np.linspace(-2, 2, 11)]}) + bbo = s.BlackBox([{'x': x} for x in np.linspace(-2, 2, 11)]) + bbo = s.BlackBox([[x] for x in np.linspace(-2, 2, 11)]) + bbo = s.BlackBox([[x] for x in np.linspace(-2, 2, 11)], True, optn=False) + bbo = s.BlackBox([[x] for x in np.linspace(-2, 2, 11)], optn1=True, optn2=False) + # END: AdvancedRTC_Snippet_01 + + def test_edi_snippet_advancedRTC_02(self): + import numpy as np + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + class SignomialTest(BlackBoxFunctionModel): + def __init__(self): + # Set up all the attributes by calling Model.__init__ + super().__init__() + + # Setup Inputs + self.inputs.append('x', '', 'Independent Variable') + + # Setup Outputs + self.outputs.append('y', '', 'Dependent Variable') + + # Simple model description + self.description = ( + 'This model evaluates the ' + 'function: max([-6*x-6, x**4-3*x**2])' + ) + + self.availableDerivative = 1 + + # BEGIN: AdvancedRTC_Snippet_02 + def BlackBox(self, *args, **kwargs): + runCases, extras = self.parseInputs(*args, **kwargs) + # END: AdvancedRTC_Snippet_02 + + x = np.array( + [pyo.value(runCases[i]['x']) for i in range(0, len(runCases))] + ) + + y = np.maximum(-6 * x - 6, x**4 - 3 * x**2) + dydx = 4 * x**3 - 6 * x + ddy_ddx = 12 * x**2 - 6 + gradientSwitch = -6 * x - 6 > x**4 - 3 * x**2 + dydx[gradientSwitch] = -6 + ddy_ddx[gradientSwitch] = 0 + + y = [yval * units.dimensionless for yval in y] + dydx = [dydx[i] * units.dimensionless for i in range(0, len(dydx))] + + returnMode = 1 + + if returnMode < 0: + returnMode = -1 * (returnMode + 1) + if returnMode == 0: + return y[0] + if returnMode == 1: + return y[0], dydx + else: + if returnMode == 0: + opt = [] + for i in range(0, len(y)): + opt.append([y[i]]) + return opt + if returnMode == 1: + opt = [] + for i in range(0, len(y)): + opt.append([[y[i]], [[[dydx[i]]]]]) + return opt + + s = SignomialTest() + ivals = [[x] for x in np.linspace(-2, 2, 11)] + + # How the black box may be called using EDI + bbo = s.BlackBox(**{'x': 0.5}) + bbo = s.BlackBox({'x': 0.5}) + bbo = s.BlackBox(**{'x': 0.5, 'optn': True}) + + # Additional options available with parseInputs + bbo = s.BlackBox(*[0.5], **{'optn1': True, 'optn2': False}) + bbo = s.BlackBox(*[0.5, True], **{'optn': False}) + bbo = s.BlackBox({'x': [x for x in np.linspace(-2, 2, 11)]}) + bbo = s.BlackBox([{'x': x} for x in np.linspace(-2, 2, 11)]) + bbo = s.BlackBox([[x] for x in np.linspace(-2, 2, 11)]) + bbo = s.BlackBox([[x] for x in np.linspace(-2, 2, 11)], True, optn=False) + bbo = s.BlackBox([[x] for x in np.linspace(-2, 2, 11)], optn1=True, optn2=False) + + +if __name__ == '__main__': + unittest.main() diff --git a/pyomo/contrib/edi/tests/test_examples.py b/pyomo/contrib/edi/tests/test_examples.py new file mode 100644 index 00000000000..72eea036838 --- /dev/null +++ b/pyomo/contrib/edi/tests/test_examples.py @@ -0,0 +1,93 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# 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. +# +# Development of this module was conducted as part of the Institute for +# the Design of Advanced Energy Systems (IDAES) with support through the +# Simulation-Based Engineering, Crosscutting Research Program within the +# U.S. Department of Energy’s Office of Fossil Energy and Carbon Management. +# +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import pyomo.common.unittest as unittest +import pyomo.environ as pyo +from pyomo.common.dependencies import attempt_import + +testIndex = 0 +import importlib + +from pyomo.core.base.units_container import pint_available + +from pyomo.common.dependencies import numpy, numpy_available +from pyomo.common.dependencies import scipy, scipy_available + +egb, egb_available = attempt_import( + "pyomo.contrib.pynumero.interfaces.external_grey_box" +) + +formulation_available = False +try: + from pyomo.contrib.edi import Formulation + + formulation_available = True +except: + pass + # formulation_available = False + +blackbox_available = False +try: + from pyomo.contrib.edi import BlackBoxFunctionModel + + blackbox_available = True +except: + pass + # blackbox_available = False + +if numpy_available: + import numpy as np + + +@unittest.skipIf( + not egb_available, 'Testing pyomo.contrib.edi requires pynumero external grey boxes' +) +@unittest.skipIf(not formulation_available, 'Formulation import failed') +@unittest.skipIf(not blackbox_available, 'Blackbox import failed') +@unittest.skipIf(not numpy_available, 'Testing pyomo.contrib.edi requires numpy') +@unittest.skipIf(not scipy_available, 'Testing pyomo.contrib.edi requires scipy') +@unittest.skipIf(not pint_available, 'Testing units requires pint') +class EDIExamples(unittest.TestCase): + def test_edi_example_placeholder(self): + "A placeholder" + pass + + +def create_new(filename): + def t_function(self): + importName = filename[0:-3] + # filename = ".."+filename + try: + importlib.import_module("pyomo.contrib.edi.examples." + importName) + except: + self.fail("This example is failing: %s" % (filename)) + + return t_function + + +pythonFileList = ["readme_example.py", "aircraft_gp.py"] + +for filename in pythonFileList: + testName = 'test_DocumentationExample_%d' % (testIndex) + testIndex += 1 + t_Function = create_new(filename) + if pint_available: + setattr(EDIExamples, testName, t_Function) + + +if __name__ == '__main__': + unittest.main() diff --git a/pyomo/contrib/edi/tests/test_formulation.py b/pyomo/contrib/edi/tests/test_formulation.py new file mode 100644 index 00000000000..37d7a9ae53e --- /dev/null +++ b/pyomo/contrib/edi/tests/test_formulation.py @@ -0,0 +1,949 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# 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. +# +# Development of this module was conducted as part of the Institute for +# the Design of Advanced Energy Systems (IDAES) with support through the +# Simulation-Based Engineering, Crosscutting Research Program within the +# U.S. Department of Energy’s Office of Fossil Energy and Carbon Management. +# +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import pyomo.common.unittest as unittest +import pyomo.environ as pyo +from pyomo.common.dependencies import attempt_import + +from pyomo.core.base.units_container import pint_available + +from pyomo.common.dependencies import numpy, numpy_available +from pyomo.common.dependencies import scipy, scipy_available + +egb, egb_available = attempt_import( + "pyomo.contrib.pynumero.interfaces.external_grey_box" +) + +formulation_available = False +try: + from pyomo.contrib.edi import Formulation + + formulation_available = True +except: + pass + # formulation_available = False + +blackbox_available = False +try: + from pyomo.contrib.edi import BlackBoxFunctionModel + + blackbox_available = True +except: + pass + # blackbox_available = False + +if numpy_available: + import numpy as np + + +@unittest.skipIf( + not egb_available, 'Testing pyomo.contrib.edi requires pynumero external grey boxes' +) +@unittest.skipIf(not formulation_available, 'Formulation import failed') +@unittest.skipIf(not blackbox_available, 'Blackbox import failed') +@unittest.skipIf(not numpy_available, 'Testing pyomo.contrib.edi requires numpy') +@unittest.skipIf(not scipy_available, 'Testing pyomo.contrib.edi requires scipy') +@unittest.skipIf(not pint_available, 'Testing units requires pint') +class TestEDIFormulation(unittest.TestCase): + def test_edi_formulation_init(self): + "Tests that a formulation initializes to the correct type and has proper data" + from pyomo.environ import ConcreteModel + from pyomo.contrib.edi import Formulation + + f = Formulation() + + self.assertIsInstance(f, Formulation) + self.assertIsInstance(f, ConcreteModel) + + self.assertEqual(f._objective_counter, 0) + self.assertEqual(f._constraint_counter, 0) + self.assertEqual(f._variable_keys, []) + self.assertEqual(f._constant_keys, []) + self.assertEqual(f._objective_keys, []) + self.assertEqual(f._runtimeObjective_keys, []) + self.assertEqual(f._objective_keys, []) + self.assertEqual(f._runtimeConstraint_keys, []) + self.assertEqual(f._constraint_keys, []) + self.assertEqual(f._allConstraint_keys, []) + + def test_edi_formulation_variable(self): + "Tests the variable constructor in edi.formulation" + import pyomo + from pyomo.contrib.edi import Formulation + from pyomo.environ import Reals, PositiveReals + + f = Formulation() + + x1 = f.Variable( + name='x1', + guess=1.0, + units='m', + description='The x variable', + size=None, + bounds=None, + domain=None, + ) + self.assertRaises(RuntimeError, f.Variable, *('x1', 1.0, 'm')) + x2 = f.Variable('x2', 1.0, 'm') + x3 = f.Variable('x3', 1.0, 'm', 'The x variable', None, None, None) + x4 = f.Variable( + name='x4', + guess=1.0, + units='m', + description='The x variable', + size=None, + bounds=None, + domain=PositiveReals, + ) + self.assertRaises( + RuntimeError, + f.Variable, + **{ + 'name': 'x5', + 'guess': 1.0, + 'units': 'm', + 'description': 'The x variable', + 'size': None, + 'bounds': None, + 'domain': "error", + } + ) + + x6 = f.Variable( + name='x6', + guess=1.0, + units='m', + description='The x variable', + size=0, + bounds=None, + domain=None, + ) + x7 = f.Variable( + name='x7', + guess=1.0, + units='m', + description='The x variable', + size=5, + bounds=None, + domain=None, + ) + self.assertRaises( + ValueError, + f.Variable, + **{ + 'name': 'x8', + 'guess': 1.0, + 'units': 'm', + 'description': 'The x variable', + 'size': 'error', + 'bounds': None, + 'domain': None, + } + ) + x9 = f.Variable( + name='x9', + guess=1.0, + units='m', + description='The x variable', + size=[2, 2], + bounds=None, + domain=None, + ) + self.assertRaises( + ValueError, + f.Variable, + **{ + 'name': 'x10', + 'guess': 1.0, + 'units': 'm', + 'description': 'The x variable', + 'size': ['2', '2'], + 'bounds': None, + 'domain': None, + } + ) + self.assertRaises( + ValueError, + f.Variable, + **{ + 'name': 'x11', + 'guess': 1.0, + 'units': 'm', + 'description': 'The x variable', + 'size': [2, 1], + 'bounds': None, + 'domain': None, + } + ) + + x12 = f.Variable( + name='x12', + guess=1.0, + units='m', + description='The x variable', + size=None, + bounds=[-10, 10], + domain=None, + ) + self.assertRaises( + ValueError, + f.Variable, + **{ + 'name': 'x13', + 'guess': 1.0, + 'units': 'm', + 'description': 'The x variable', + 'size': None, + 'bounds': [10, -10], + 'domain': None, + } + ) + self.assertRaises( + ValueError, + f.Variable, + **{ + 'name': 'x14', + 'guess': 1.0, + 'units': 'm', + 'description': 'The x variable', + 'size': None, + 'bounds': ["-10", "10"], + 'domain': None, + } + ) + self.assertRaises( + ValueError, + f.Variable, + **{ + 'name': 'x15', + 'guess': 1.0, + 'units': 'm', + 'description': 'The x variable', + 'size': None, + 'bounds': [1, 2, 3], + 'domain': None, + } + ) + self.assertRaises( + ValueError, + f.Variable, + **{ + 'name': 'x16', + 'guess': 1.0, + 'units': 'm', + 'description': 'The x variable', + 'size': None, + 'bounds': "error", + 'domain': None, + } + ) + self.assertRaises( + ValueError, + f.Variable, + **{ + 'name': 'x17', + 'guess': 1.0, + 'units': 'm', + 'description': 'The x variable', + 'size': None, + 'bounds': [0, "10"], + 'domain': None, + } + ) + + # verifies alternate unit construction + x18 = f.Variable('x18', 1.0, pyo.units.m) + self.assertRaises(AttributeError, f.Variable, *('x19', 1.0, 'string')) + + def test_edi_formulation_constant(self): + "Tests the constant constructor in edi.formulation" + from pyomo.contrib.edi import Formulation + from pyomo.environ import Reals, PositiveReals + + f = Formulation() + + c1 = f.Constant( + name='c1', + value=1.0, + units='m', + description='A constant c', + size=None, + within=None, + ) + self.assertRaises(RuntimeError, f.Constant, *('c1', 1.0, 'm')) + c2 = f.Constant('c2', 1.0, 'm') + c3 = f.Constant('c3', 1.0, 'm', 'A constant c', None, None) + c4 = f.Constant( + name='c4', + value=1.0, + units='m', + description='A constant c', + size=None, + within=PositiveReals, + ) + self.assertRaises( + RuntimeError, + f.Constant, + **{ + 'name': 'c5', + 'value': 1.0, + 'units': 'm', + 'description': 'The x variable', + 'size': None, + 'within': "error", + } + ) + + c6 = f.Constant( + name='c6', + value=1.0, + units='m', + description='A constant c', + size=0, + within=None, + ) + c7 = f.Constant( + name='c7', + value=1.0, + units='m', + description='A constant c', + size=5, + within=None, + ) + self.assertRaises( + ValueError, + f.Constant, + **{ + 'name': 'c8', + 'value': 1.0, + 'units': 'm', + 'description': 'A constant c', + 'size': 'error', + 'within': None, + } + ) + c9 = f.Constant( + name='c9', + value=1.0, + units='m', + description='A constant c', + size=[2, 2], + within=None, + ) + self.assertRaises( + ValueError, + f.Constant, + **{ + 'name': 'c10', + 'value': 1.0, + 'units': 'm', + 'description': 'A constant c', + 'size': ['2', '2'], + 'within': None, + } + ) + self.assertRaises( + ValueError, + f.Constant, + **{ + 'name': 'c11', + 'value': 1.0, + 'units': 'm', + 'description': 'A constant c', + 'size': [2, 1], + 'within': None, + } + ) + + def test_edi_formulation_objective(self): + "Tests the objective constructor in edi.formulation" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + f.Objective(x + y) + + def test_edi_formulation_runtimeobjective(self): + "Tests the runtime objective constructor in edi.formulation" + # TODO: not currently implemented, see: https://github.com/codykarcher/pyomo/issues/5 + pass + + def test_edi_formulation_constraint(self): + "Tests the constraint constructor in edi.formulation" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + f.Objective(x + y) + f.Constraint(x + y <= 1.0 * units.m) + + def test_edi_formulation_runtimeconstraint_tuple(self): + "Tests the runtime constraint constructor in edi.formulation" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + z = f.Variable( + name='z', guess=1.0, units='m^2', description='The unit circle output' + ) + c = f.Constant( + name='c', value=1.0, units='', description='A constant c', size=2 + ) + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super(UnitCircle, self).__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', + units='ft**2', + description='Resultant of the unit circle evaluation', + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + x = pyo.value( + units.convert(x, self.inputs[0].units) + ) # Converts to correct units then casts to float + y = pyo.value( + units.convert(y, self.inputs[1].units) + ) # Converts to correct units then casts to float + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + f.Constraint(z <= 1 * units.m**2) + + f.RuntimeConstraint(*(z, '==', [x, y], UnitCircle())) + + def test_edi_formulation_runtimeconstraint_list(self): + "Tests the runtime constraint constructor in edi.formulation" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + z = f.Variable( + name='z', guess=1.0, units='m^2', description='The unit circle output' + ) + c = f.Constant( + name='c', value=1.0, units='', description='A constant c', size=2 + ) + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super(UnitCircle, self).__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', + units='ft**2', + description='Resultant of the unit circle evaluation', + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + x = pyo.value( + units.convert(x, self.inputs[0].units) + ) # Converts to correct units then casts to float + y = pyo.value( + units.convert(y, self.inputs[1].units) + ) # Converts to correct units then casts to float + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + f.Constraint(z <= 1 * units.m**2) + + f.RuntimeConstraint(*[[z], ['=='], [x, y], UnitCircle()]) + + def test_edi_formulation_runtimeconstraint_dict(self): + "Tests the runtime constraint constructor in edi.formulation" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + z = f.Variable( + name='z', guess=1.0, units='m^2', description='The unit circle output' + ) + c = f.Constant( + name='c', value=1.0, units='', description='A constant c', size=2 + ) + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super(UnitCircle, self).__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', + units='ft**2', + description='Resultant of the unit circle evaluation', + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + x = pyo.value( + units.convert(x, self.inputs[0].units) + ) # Converts to correct units then casts to float + y = pyo.value( + units.convert(y, self.inputs[1].units) + ) # Converts to correct units then casts to float + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + f.Constraint(z <= 1 * units.m**2) + + f.RuntimeConstraint( + **{ + 'outputs': z, + 'operators': '==', + 'inputs': [x, y], + 'black_box': UnitCircle(), + } + ) + + def test_edi_formulation_constraintlist_1(self): + "Tests the constraint list constructor in edi.formulation" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + z = f.Variable( + name='z', guess=1.0, units='m^2', description='The unit circle output' + ) + c = f.Constant( + name='c', value=1.0, units='', description='A constant c', size=2 + ) + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super(UnitCircle, self).__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', + units='ft**2', + description='Resultant of the unit circle evaluation', + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + x = pyo.value( + units.convert(x, self.inputs[0].units) + ) # Converts to correct units then casts to float + y = pyo.value( + units.convert(y, self.inputs[1].units) + ) # Converts to correct units then casts to float + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + f.ConstraintList([(z, '==', [x, y], UnitCircle()), z <= 1 * units.m**2]) + + cl = f.get_constraints() + + self.assertTrue(len(cl) == 2) + + def test_edi_formulation_constraintlist_2(self): + "Tests the constraint list constructor in edi.formulation" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m**2', description='The y variable') + f.Objective(y) + + class Parabola(BlackBoxFunctionModel): + def __init__(self): + super(Parabola, self).__init__() + self.description = 'This model evaluates the function: y = x**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.outputs.append( + name='y', units='ft**2', description='The y variable' + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x): # The actual function that does things + x = pyo.value( + units.convert(x, self.inputs[0].units) + ) # Converts to correct units then casts to float + y = x**2 # Compute y + dydx = 2 * x # Compute dy/dx + y *= units.ft**2 + dydx *= units.ft # units.ft**2 / units.ft + return y, [dydx] # return z, grad(z), hess(z)... + + f.ConstraintList( + [{'outputs': y, 'operators': '==', 'inputs': x, 'black_box': Parabola()}] + ) + + cl = f.get_constraints() + + self.assertTrue(len(cl) == 1) + + def test_edi_formulation_constraintlist_3(self): + "Tests the constraint list constructor in edi.formulation" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable( + name='x', guess=1.0, units='m', description='The x variable', size=3 + ) + y = f.Variable( + name='y', guess=1.0, units='m**2', description='The y variable', size=3 + ) + f.Objective(y[0] + y[1] + y[2]) + + class Parabola(BlackBoxFunctionModel): + def __init__(self): + super(Parabola, self).__init__() + self.description = 'This model evaluates the function: y = x**2' + self.inputs.append( + name='x', size=3, units='ft', description='The x variable' + ) + self.outputs.append( + name='y', size=3, units='ft**2', description='The y variable' + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(*args, **kwargs): # The actual function that does things + args = list(args) + self = args.pop(0) + runCases, returnMode, remainingKwargs = self.parseInputs( + *args, **kwargs + ) + + x = self.sanitizeInputs(runCases[0]['x']) + x = np.array([pyo.value(xval) for xval in x], dtype=np.float64) + + y = x**2 # Compute y + dydx = 2 * x # Compute dy/dx + + y = y * units.ft**2 + dydx = np.diag(dydx) + dydx = dydx * units.ft # units.ft**2 / units.ft + + return y, [dydx] # return z, grad(z), hess(z)... + + f.ConstraintList( + [{'outputs': y, 'operators': '==', 'inputs': x, 'black_box': Parabola()}] + ) + cl = f.get_constraints() + self.assertTrue(len(cl) == 1) + + def test_edi_formulation_runtimeconstraint_exceptions(self): + "Tests the runtime constraint constructor in edi.formulation" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + z = f.Variable( + name='z', guess=1.0, units='m^2', description='The unit circle output' + ) + c = f.Constant( + name='c', value=1.0, units='', description='A constant c', size=2 + ) + f.Objective(x + y) + + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): + super(UnitCircle, self).__init__() + self.description = 'This model evaluates the function: z = x**2 + y**2' + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + self.outputs.append( + name='z', + units='ft**2', + description='Resultant of the unit circle evaluation', + ) + self.availableDerivative = 1 + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + x = pyo.value( + units.convert(x, self.inputs[0].units) + ) # Converts to correct units then casts to float + y = pyo.value( + units.convert(y, self.inputs[1].units) + ) # Converts to correct units then casts to float + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + f.Constraint(z <= 1 * units.m**2) + + # flaggs the input or as bad before assigning the incorrect black box + self.assertRaises( + ValueError, f.RuntimeConstraint, *(z, '==', 1.0, UnitCircle()) + ) + self.assertRaises( + ValueError, f.RuntimeConstraint, *(1.0, '==', [x, y], UnitCircle()) + ) + + self.assertRaises( + ValueError, f.RuntimeConstraint, *(z, '==', [1.0, y], UnitCircle()) + ) + self.assertRaises( + ValueError, f.RuntimeConstraint, *(z, 1.0, [x, y], UnitCircle()) + ) + self.assertRaises( + ValueError, f.RuntimeConstraint, *(z, '=', [x, y], UnitCircle()) + ) + + def test_edi_formulation_getvariables(self): + "Tests the get_variables function in edi.formulation" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + + vrs = f.get_variables() + self.assertListEqual(vrs, [x, y]) + + def test_edi_formulation_getconstants(self): + "Tests the get_constants function in edi.formulation" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + + c1 = f.Constant( + name='c1', + value=1.0, + units='m', + description='A constant c1', + size=None, + within=None, + ) + c2 = f.Constant( + name='c2', + value=1.0, + units='m', + description='A constant c2', + size=None, + within=None, + ) + + csts = f.get_constants() + self.assertListEqual(csts, [c1, c2]) + + def test_edi_formulation_getobjectives(self): + "Tests the get_objectives function in edi.formulation" + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + f.Objective(x + y) + objList = f.get_objectives() + self.assertTrue(len(objList) == 1) + # not really sure how to check this, so I wont + + def test_edi_formulation_getconstraints(self): + "Tests the get_constraints, get_explicitConstraints, and get_runtimeConstraints functions in edi.formulation" + # ================= + # Import Statements + # ================= + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation, BlackBoxFunctionModel + + # =================== + # Declare Formulation + # =================== + f = Formulation() + + # ================= + # Declare Variables + # ================= + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + z = f.Variable( + name='z', guess=1.0, units='m^2', description='The unit circle output' + ) + + # ================= + # Declare Constants + # ================= + c = f.Constant( + name='c', value=1.0, units='', description='A constant c', size=2 + ) + + # ===================== + # Declare the Objective + # ===================== + f.Objective(c[0] * x + c[1] * y) + + # =================== + # Declare a Black Box + # =================== + class UnitCircle(BlackBoxFunctionModel): + def __init__(self): # The initialization function + # Initialize the black box model + super(UnitCircle, self).__init__() + + # A brief description of the model + self.description = 'This model evaluates the function: z = x**2 + y**2' + + # Declare the black box model inputs + self.inputs.append(name='x', units='ft', description='The x variable') + self.inputs.append(name='y', units='ft', description='The y variable') + + # Declare the black box model outputs + self.outputs.append( + name='z', + units='ft**2', + description='Resultant of the unit circle evaluation', + ) + + # Declare the maximum available derivative + self.availableDerivative = 1 + + # Post-initialization setup + self.post_init_setup(len(self.inputs)) + + def BlackBox(self, x, y): # The actual function that does things + x = pyo.value( + units.convert(x, self.inputs[0].units) + ) # Converts to correct units then casts to float + y = pyo.value( + units.convert(y, self.inputs[1].units) + ) # Converts to correct units then casts to float + + z = x**2 + y**2 # Compute z + dzdx = 2 * x # Compute dz/dx + dzdy = 2 * y # Compute dz/dy + + z *= units.ft**2 + dzdx *= units.ft # units.ft**2 / units.ft + dzdy *= units.ft # units.ft**2 / units.ft + + return z, [dzdx, dzdy] # return z, grad(z), hess(z)... + + # ======================= + # Declare the Constraints + # ======================= + f.ConstraintList([(z, '==', [x, y], UnitCircle()), z <= 1 * units.m**2]) + + cl = f.get_constraints() + ecl = f.get_explicitConstraints() + rcl = f.get_runtimeConstraints() + + self.assertTrue(len(cl) == 2) + self.assertTrue(len(ecl) == 1) + self.assertTrue(len(rcl) == 1) + + def test_edi_formulation_checkunits(self): + "Tests the check_units function in edi.formulation" + import pyomo + import pyomo.environ as pyo + from pyomo.environ import units + from pyomo.contrib.edi import Formulation + + f = Formulation() + x = f.Variable(name='x', guess=1.0, units='m', description='The x variable') + y = f.Variable(name='y', guess=1.0, units='m', description='The y variable') + + f.Objective(x + y) + f.Constraint(x + y <= 1.0 * units.m) + f.check_units() + + f.Constraint(2.0 * x + y <= 1.0) + self.assertRaises( + pyomo.core.base.units_container.UnitsError, f.check_units, *() + ) + + f2 = Formulation() + u = f2.Variable(name='u', guess=1.0, units='m', description='The u variable') + v = f2.Variable(name='v', guess=1.0, units='kg', description='The v variable') + f2.Objective(u + v) + self.assertRaises( + pyomo.core.base.units_container.UnitsError, f2.check_units, *() + ) + + +if __name__ == '__main__': + unittest.main()