pyomo

---
name: pyomo
description: Expert in pyomo modelling in python
---

# Pyomo Skill

Comprehensive assistance with pyomo development, generated from official documentation.

## When to Use This Skill

This skill should be triggered when:
- Working with pyomo
- Asking about pyomo features or APIs
- Implementing pyomo solutions
- Debugging pyomo code
- Learning pyomo best practices

## Quick Reference

### Common Patterns

**Pattern 1:** One can add terms to an objective function of a ConcreteModel (or and instantiated AbstractModel) using the expr attribute of the objective function object. Here is a simple example:

```
ConcreteModel
```

**Pattern 2:** Multiple objectives can be declared, but only one can be active at a time (at present, Pyomo does not support any solvers that can be given more than one objective). If both model.obj1 and model.obj2 have been declared using Objective, then one can ensure that model.obj2 is passed to the solver as shown in this simple example:

```
model.obj1
```

**Pattern 3:** The Pyomo Configuration System The Pyomo configuration system provides a set of three classes (ConfigDict, ConfigList, and ConfigValue) for managing and documenting structured configuration information and user input. The system is based around the ConfigValue class, which provides storage for a single configuration entry. ConfigValue objects can be grouped using two containers (ConfigDict and ConfigList), which provide functionality analogous to Python’s dict and list classes, respectively. At its simplest, the configuration system allows for developers to specify a dictionary of documented configuration entries: from pyomo.common.config import ( ConfigDict, ConfigList, ConfigValue ) config = ConfigDict() config.declare('filename', ConfigValue( default=None, domain=str, description="Input file name", )) config.declare("bound tolerance", ConfigValue( default=1E-5, domain=float, description="Bound tolerance", doc="Relative tolerance for bound feasibility checks" )) config.declare("iteration limit", ConfigValue( default=30, domain=int, description="Iteration limit", doc="Number of maximum iterations in the decomposition methods" )) Users can then provide values for those entries, and retrieve the current values: >>> config['filename'] = 'tmp.txt' >>> print(config['filename']) tmp.txt >>> print(config['iteration limit']) 30 For convenience, ConfigDict objects support read/write access via attributes (with spaces in the declaration names replaced by underscores): >>> print(config.filename) tmp.txt >>> print(config.iteration_limit) 30 >>> config.iteration_limit = 20 >>> print(config.iteration_limit) 20 Domain validation All Config objects support a domain keyword that accepts a callable object (type, function, or callable instance). The domain callable should take data and map it onto the desired domain, optionally performing domain validation (see ConfigValue, ConfigDict, and ConfigList for more information). This allows client code to accept a very flexible set of inputs without “cluttering” the code with input validation: >>> config.iteration_limit = 35.5 >>> print(config.iteration_limit) 35 >>> print(type(config.iteration_limit).__name__) int In addition to common types (like int, float, bool, and str), the configuration system provides a number of custom domain validators for common use cases: Bool(val) Domain validator for bool-like objects. Integer(val) Domain validation function admitting integers PositiveInt(val) Domain validation function admitting strictly positive integers NegativeInt(val) Domain validation function admitting strictly negative integers NonNegativeInt(val) Domain validation function admitting integers >= 0 NonPositiveInt(val) Domain validation function admitting integers <= 0 PositiveFloat(val) Domain validation function admitting strictly positive numbers NegativeFloat(val) Domain validation function admitting strictly negative numbers NonPositiveFloat(val) Domain validation function admitting numbers less than or equal to 0 NonNegativeFloat(val) Domain validation function admitting numbers greater than or equal to 0 In(domain[, cast]) Domain validation class admitting a Container of possible values InEnum(domain) Domain validation class admitting an enum value/name. IsInstance(*bases[, document_full_base_names]) Domain validator for type checking. ListOf(itemtype[, domain, string_lexer]) Domain validator for lists of a specified type Module([basePath, expandPath]) Domain validator for modules. Path([basePath, expandPath]) Domain validator for a path-like object. PathList([basePath, expandPath]) Domain validator for a list of path-like objects. DynamicImplicitDomain(callback) Implicit domain that can return a custom domain based on the key. Configuring class hierarchies A feature of the configuration system is that the core classes all implement __call__, and can themselves be used as domain values. Beyond providing domain verification for complex hierarchical structures, this feature allows ConfigDict objects to cleanly support extension and the configuration of derived classes. Consider the following example: >>> class Base: ... CONFIG = ConfigDict() ... CONFIG.declare('filename', ConfigValue( ... default='input.txt', ... domain=str, ... )) ... def __init__(self, **kwds): ... self.cfg = self.CONFIG(kwds) ... self.cfg.display() ... >>> class Derived(Base): ... CONFIG = Base.CONFIG() ... CONFIG.declare('pattern', ConfigValue( ... default=None, ... domain=str, ... )) ... >>> tmp = Base(filename='foo.txt') filename: foo.txt >>> tmp = Derived(pattern='.*warning') filename: input.txt pattern: .*warning Here, the base class Base declares a class-level attribute CONFIG as a ConfigDict containing a single entry (filename). The derived class (Derived) then starts by making a copy of the base class’ CONFIG, and then defines an additional entry (pattern). Instances of the base class will still create cfg attributes that only have the single filename entry, whereas instances of the derived class will have cfg attributes with two entries: the pattern entry declared by the derived class, and the filename entry “inherited” from the base class. An extension of this design pattern provides a clean approach for handling “ephemeral” instance options. Consider an interface to an external “solver”. Our class implements a solve() method that takes a problem and sends it to the solver along with some solver configuration options. We would like to be able to set those options “persistently” on instances of the interface class, but still override them “temporarily” for individual calls to solve(). We implement this by creating copies of the class’s configuration for both specific instances and for use by each solve() call: class Solver: CONFIG = ConfigDict() CONFIG.declare('iterlim', ConfigValue( default=10, domain=int, )) def __init__(self, **kwds): self.config = self.CONFIG(kwds) def solve(self, model, **options): config = self.config(options) # Solve the model with the specified iterlim config.display() >>> solver = Solver() >>> solver.solve(None) iterlim: 10 >>> solver.config.iterlim = 20 >>> solver.solve(None) iterlim: 20 >>> solver.solve(None, iterlim=50) iterlim: 50 >>> solver.solve(None) iterlim: 20 This design pattern is widely used across Pyomo; particularly for configuring solver interfaces and transformations. We provide a decorator to simplify the process of documenting these CONFIG attributes: from pyomo.common.config import document_class_CONFIG @document_class_CONFIG(methods=['solve']) class MySolver: """Interface to My Solver""" # #: Global class configuration; see :ref:`MySolver_CONFIG` CONFIG = ConfigDict() CONFIG.declare('iterlim', ConfigValue( default=10, domain=int, doc="Solver iteration limit", )) # def __init__(self, **kwds): #: Instance configuration; see :ref:`MySolver_CONFIG` self.config = self.CONFIG(kwds) # def solve(self, model, **options): """Solve `model` using My Solver""" # config = self.config(options) # Solve the model with the specified iterlim config.display() >>> print(MySolver.__doc__) Interface to My Solver **Class configuration** This class leverages the Pyomo Configuration System for managing configuration options. See the discussion on :ref:`configuring class hierarchies <class_config>` for more information on how configuration class attributes, instance attributes, and method keyword arguments interact. .. _MySolver::CONFIG: CONFIG ------ iterlim: int, default=10 Solver iteration limit >>> print(MySolver.solve.__doc__) Solve `model` using My Solver Keyword Arguments ----------------- iterlim: int, default=10 Solver iteration limit Interacting with argparse In addition to basic storage and retrieval, the configuration system provides hooks to the argparse command-line argument parsing system. Individual configuration entries can be declared as argparse arguments using the declare_as_argument() method. To make declaration simpler, the declare() method returns the declared configuration object so that the argument declaration can be done inline: import argparse config = ConfigDict() config.declare('iterlim', ConfigValue( domain=int, default=100, description="iteration limit", )).declare_as_argument() config.declare('lbfgs', ConfigValue( domain=bool, description="use limited memory BFGS update", )).declare_as_argument() config.declare('linesearch', ConfigValue( domain=bool, default=True, description="use line search", )).declare_as_argument() config.declare('relative tolerance', ConfigValue( domain=float, description="relative convergence tolerance", )).declare_as_argument('--reltol', '-r', group='Tolerances') config.declare('absolute tolerance', ConfigValue( domain=float, description="absolute convergence tolerance", )).declare_as_argument('--abstol', '-a', group='Tolerances') The ConfigDict can then be used to initialize (or augment) an argparse.ArgumentParser object: parser = argparse.ArgumentParser("tester") config.initialize_argparse(parser) Key information from the ConfigDict is automatically transferred over to the ArgumentParser object: >>> print(parser.format_help()) usage: tester [-h] [--iterlim INT] [--lbfgs] [--disable-linesearch] [--reltol FLOAT] [--abstol FLOAT] ... -h, --help show this help message and exit --iterlim INT iteration limit --lbfgs use limited memory BFGS update --disable-linesearch [DON'T] use line search Tolerances: --reltol... -r FLOAT relative convergence tolerance --abstol... -a FLOAT absolute convergence tolerance Parsed arguments can then be imported back into the ConfigDict: >>> args=parser.parse_args(['--lbfgs', '--reltol', '0.1', '-a', '0.2']) >>> args = config.import_argparse(args) >>> config.display() iterlim: 100 lbfgs: true linesearch: true relative tolerance: 0.1 absolute tolerance: 0.2 Accessing user-specified values It is frequently useful to know which values a user explicitly set, and which values a user explicitly set but have never been retrieved. The configuration system provides two generator methods to return the items that a user explicitly set (user_values()) and the items that were set but never retrieved (unused_user_values()): >>> print([val.name() for val in config.user_values()]) ['lbfgs', 'relative tolerance', 'absolute tolerance'] >>> print(config.relative_tolerance) 0.1 >>> print([val.name() for val in config.unused_user_values()]) ['lbfgs', 'absolute tolerance'] Generating output & documentation Configuration objects support three methods for generating output and documentation: display(), generate_yaml_template(), and generate_documentation(). The simplest is display(), which prints out the current values of the configuration object (and if it is a container type, all of its children). generate_yaml_template() is similar to display(), but also includes the description fields as formatted comments. solver_config = config config = ConfigDict() config.declare('output', ConfigValue( default='results.yml', domain=str, description='output results filename' )) config.declare('verbose', ConfigValue( default=0, domain=int, description='output verbosity', doc='This sets the system verbosity. The default (0) only logs ' 'warnings and errors. Larger integer values will produce ' 'additional log messages.', )) config.declare('solvers', ConfigList( domain=solver_config, description='list of solvers to apply', )) >>> config.display() output: results.yml verbose: 0 solvers: [] >>> print(config.generate_yaml_template()) output: results.yml # output results filename verbose: 0 # output verbosity solvers: [] # list of solvers to apply It is important to note that both methods document the current state of the configuration object. So, in the example above, since the solvers list is empty, you will not get any information on the elements in the list. Of course, if you add a value to the list, then the data will be output: >>> tmp = config() >>> tmp.solvers.append({}) >>> tmp.display() output: results.yml verbose: 0 solvers: - iterlim: 100 lbfgs: true linesearch: true relative tolerance: 0.1 absolute tolerance: 0.2 >>> print(tmp.generate_yaml_template()) output: results.yml # output results filename verbose: 0 # output verbosity solvers: # list of solvers to apply - iterlim: 100 # iteration limit lbfgs: true # use limited memory BFGS update linesearch: true # use line search relative tolerance: 0.1 # relative convergence tolerance absolute tolerance: 0.2 # absolute convergence tolerance The third method (generate_documentation()) behaves differently. This method is designed to generate reference documentation. For each configuration item, the doc field is output. If the item has no doc, then the description field is used. List containers have their domain documented and not their current values. The documentation can be configured through optional arguments. The defaults generate LaTeX documentation: >>> print(config.generate_documentation()) \begin{description}[topsep=0pt,parsep=0.5em,itemsep=-0.4em] \item[{output}]\hfill \\output results filename \item[{verbose}]\hfill \\This sets the system verbosity. The default (0) only logs warnings and errors. Larger integer values will produce additional log messages. \item[{solvers}]\hfill \\list of solvers to apply \begin{description}[topsep=0pt,parsep=0.5em,itemsep=-0.4em] \item[{iterlim}]\hfill \\iteration limit \item[{lbfgs}]\hfill \\use limited memory BFGS update \item[{linesearch}]\hfill \\use line search \item[{relative tolerance}]\hfill \\relative convergence tolerance \item[{absolute tolerance}]\hfill \\absolute convergence tolerance \end{description} \end{description}

```
ConfigDict
```

**Pattern 4:** A feature of the configuration system is that the core classes all implement __call__, and can themselves be used as domain values. Beyond providing domain verification for complex hierarchical structures, this feature allows ConfigDict objects to cleanly support extension and the configuration of derived classes. Consider the following example:

```
__call__
```

**Pattern 5:** Dynamic Optimization with pyomo.DAE The pyomo.DAE modeling extension [PyomoDAE-paper] allows users to incorporate systems of differential algebraic equations (DAE)s in a Pyomo model. The modeling components in this extension are able to represent ordinary or partial differential equations. The differential equations do not have to be written in a particular format and the components are flexible enough to represent higher-order derivatives or mixed partial derivatives. Pyomo.DAE also includes model transformations which use simultaneous discretization approaches to transform a DAE model into an algebraic model. Finally, pyomo.DAE includes utilities for simulating DAE models and initializing dynamic optimization problems. Modeling Components Pyomo.DAE introduces three new modeling components to Pyomo: pyomo.dae.ContinuousSet Represents a bounded continuous domain pyomo.dae.DerivativeVar Represents derivatives in a model and defines how a Var is differentiated pyomo.dae.Integral Represents an integral over a continuous domain As will be shown later, differential equations can be declared using using these new modeling components along with the standard Pyomo Var and Constraint components. ContinuousSet This component is used to define continuous bounded domains (for example ‘spatial’ or ‘time’ domains). It is similar to a Pyomo Set component and can be used to index things like variables and constraints. Any number of ContinuousSets can be used to index a component and components can be indexed by both Sets and ContinuousSets in arbitrary order. In the current implementation, models with ContinuousSet components may not be solved until every ContinuousSet has been discretized. Minimally, a ContinuousSet must be initialized with two numeric values representing the upper and lower bounds of the continuous domain. A user may also specify additional points in the domain to be used as finite element points in the discretization. class pyomo.dae.ContinuousSet(*args, **kwds)[source] Represents a bounded continuous domain Minimally, this set must contain two numeric values defining the bounds of a continuous range. Discrete points of interest may be added to the continuous set. A continuous set is one dimensional and may only contain numerical values. Parameters: initialize (list) – Default discretization points to be included bounds (tuple) – The bounding points for the continuous domain. The bounds will be included as discrete points in the ContinuousSet and will be used to bound the points added to the ContinuousSet through the ‘initialize’ argument, a data file, or the add() method _changed This keeps track of whether or not the ContinuousSet was changed during discretization. If the user specifies all of the needed discretization points before the discretization then there is no need to go back through the model and reconstruct things indexed by the ContinuousSet Type: boolean _fe This is a sorted list of the finite element points in the ContinuousSet. i.e. this list contains all the discrete points in the ContinuousSet that are not collocation points. Points that are both finite element points and collocation points will be included in this list. Type: list _discretization_info This is a dictionary which contains information on the discretization transformation which has been applied to the ContinuousSet. Type: dict construct(values=None)[source] Constructs a ContinuousSet component find_nearest_index(target, tolerance=None)[source] Returns the index of the nearest point in the ContinuousSet. If a tolerance is specified, the index will only be returned if the distance between the target and the closest point is less than or equal to that tolerance. If there is a tie for closest point, the index on the left is returned. Parameters: target (float) tolerance (float or None) Return type: float or None get_changed()[source] Returns flag indicating if the ContinuousSet was changed during discretization Returns “True” if additional points were added to the ContinuousSet while applying a discretization scheme Return type: boolean get_discretization_info()[source] Returns a dict with information on the discretization scheme that has been applied to the ContinuousSet. Return type: dict get_finite_elements()[source] Returns the finite element points If the ContinuousSet has been discretizaed using a collocation scheme, this method will return a list of the finite element discretization points but not the collocation points within each finite element. If the ContinuousSet has not been discretized or a finite difference discretization was used, this method returns a list of all the discretization points in the ContinuousSet. Return type: list of floats get_lower_element_boundary(point)[source] Returns the first finite element point that is less than or equal to ‘point’ Parameters: point (float) Return type: float get_upper_element_boundary(point)[source] Returns the first finite element point that is greater or equal to ‘point’ Parameters: point (float) Return type: float set_changed(newvalue)[source] Sets the _changed flag to ‘newvalue’ Parameters: newvalue (boolean) The following code snippet shows examples of declaring a ContinuousSet component on a concrete Pyomo model: Required imports >>> import pyomo.environ as pyo >>> from pyomo.dae import ContinuousSet >>> model = pyo.ConcreteModel() Declaration by providing bounds >>> model.t = ContinuousSet(bounds=(0,5)) Declaration by initializing with desired discretization points >>> model.x = ContinuousSet(initialize=[0,1,2,5]) Note A ContinuousSet may not be constructed unless at least two numeric points are provided to bound the continuous domain. The following code snippet shows an example of declaring a ContinuousSet component on an abstract Pyomo model using the example data file. set t := 0 0.5 2.25 3.75 5; Required imports >>> import pyomo.environ as pyo >>> from pyomo.dae import ContinuousSet >>> model = pyo.AbstractModel() The ContinuousSet below will be initialized using the points in the data file when a model instance is created. >>> model.t = ContinuousSet() Note If a separate data file is used to initialize a ContinuousSet, it is done using the ‘set’ command and not ‘continuousset’ Note Most valid ways to declare and initialize a Set can be used to declare and initialize a ContinuousSet. See the documentation for Set for additional options. Warning Be careful using a ContinuousSet as an implicit index in an expression, i.e. sum(m.v[i] for i in m.myContinuousSet). The expression will be generated using the discretization points contained in the ContinuousSet at the time the expression was constructed and will not be updated if additional points are added to the set during discretization. Note ContinuousSet components are always ordered (sorted) therefore the first() and last() Set methods can be used to access the lower and upper boundaries of the ContinuousSet respectively DerivativeVar class pyomo.dae.DerivativeVar(*args, **kwargs)[source] Represents derivatives in a model and defines how a Var is differentiated The DerivativeVar component is used to declare a derivative of a Var. The constructor accepts a single positional argument which is the Var that’s being differentiated. A Var may only be differentiated with respect to a ContinuousSet that it is indexed by. The indexing sets of a DerivativeVar are identical to those of the Var it is differentiating. Parameters: sVar (pyomo.environ.Var) – The variable being differentiated wrt (pyomo.dae.ContinuousSet or tuple) – Equivalent to withrespectto keyword argument. The ContinuousSet that the derivative is being taken with respect to. Higher order derivatives are represented by including the ContinuousSet multiple times in the tuple sent to this keyword. i.e. wrt=(m.t, m.t) would be the second order derivative with respect to m.t get_continuousset_list()[source] Return the a list of ContinuousSet components the derivative is being taken with respect to. Return type: list get_derivative_expression()[source] Returns the current discretization expression for this derivative or creates an access function to its Var the first time this method is called. The expression gets built up as the discretization transformations are sequentially applied to each ContinuousSet in the model. get_state_var()[source] Return the Var that is being differentiated. Return type: Var is_fully_discretized()[source] Check to see if all the ContinuousSets this derivative is taken with respect to have been discretized. Return type: boolean set_derivative_expression(expr)[source] Sets``_expr``, an expression representing the discretization equations linking the DerivativeVar to its state Var The code snippet below shows examples of declaring DerivativeVar components on a Pyomo model. In each case, the variable being differentiated is supplied as the only positional argument and the type of derivative is specified using the ‘wrt’ (or the more verbose ‘withrespectto’) keyword argument. Any keyword argument that is valid for a Pyomo Var component may also be specified. Required imports >>> import pyomo.environ as pyo >>> from pyomo.dae import ContinuousSet, DerivativeVar >>> model = pyo.ConcreteModel() >>> model.s = pyo.Set(initialize=['a','b']) >>> model.t = ContinuousSet(bounds=(0,5)) >>> model.l = ContinuousSet(bounds=(-10,10)) >>> model.x = pyo.Var(model.t) >>> model.y = pyo.Var(model.s,model.t) >>> model.z = pyo.Var(model.t,model.l) Declare the first derivative of model.x with respect to model.t >>> model.dxdt = DerivativeVar(model.x, withrespectto=model.t) Declare the second derivative of model.y with respect to model.t Note that this DerivativeVar will be indexed by both model.s and model.t >>> model.dydt2 = DerivativeVar(model.y, wrt=(model.t,model.t)) Declare the partial derivative of model.z with respect to model.l Note that this DerivativeVar will be indexed by both model.t and model.l >>> model.dzdl = DerivativeVar(model.z, wrt=(model.l), initialize=0) Declare the mixed second order partial derivative of model.z with respect to model.t and model.l and set bounds >>> model.dz2 = DerivativeVar(model.z, wrt=(model.t, model.l), bounds=(-10, 10)) Note The ‘initialize’ keyword argument will initialize the value of a derivative and is not the same as specifying an initial condition. Initial or boundary conditions should be specified using a Constraint or ConstraintList or by fixing the value of a Var at a boundary point. Declaring Differential Equations A differential equations is declared as a standard Pyomo Constraint and is not required to have any particular form. The following code snippet shows how one might declare an ordinary or partial differential equation. Required imports >>> import pyomo.environ as pyo >>> from pyomo.dae import ContinuousSet, DerivativeVar, Integral >>> model = pyo.ConcreteModel() >>> model.s = pyo.Set(initialize=['a', 'b']) >>> model.t = ContinuousSet(bounds=(0, 5)) >>> model.l = ContinuousSet(bounds=(-10, 10)) >>> model.x = pyo.Var(model.s, model.t) >>> model.y = pyo.Var(model.t, model.l) >>> model.dxdt = DerivativeVar(model.x, wrt=model.t) >>> model.dydt = DerivativeVar(model.y, wrt=model.t) >>> model.dydl2 = DerivativeVar(model.y, wrt=(model.l, model.l)) An ordinary differential equation >>> def _ode_rule(m, s, t): ... if t == 0: ... return pyo.Constraint.Skip ... return m.dxdt[s, t] == m.x[s, t]**2 >>> model.ode = pyo.Constraint(model.s, model.t, rule=_ode_rule) A partial differential equation >>> def _pde_rule(m, t, l): ... if t == 0 or l == m.l.first() or l == m.l.last(): ... return pyo.Constraint.Skip ... return m.dydt[t, l] == m.dydl2[t, l] >>> model.pde = pyo.Constraint(model.t, model.l, rule=_pde_rule) By default, a Constraint declared over a ContinuousSet will be applied at every discretization point contained in the set. Often a modeler does not want to enforce a differential equation at one or both boundaries of a continuous domain. This may be addressed explicitly in the Constraint declaration using Constraint.Skip as shown above. Alternatively, the desired constraints can be deactivated just before the model is sent to a solver as shown below. >>> def _ode_rule(m, s, t): ... return m.dxdt[s, t] == m.x[s, t]**2 >>> model.ode = pyo.Constraint(model.s, model.t, rule=_ode_rule) >>> def _pde_rule(m, t, l): ... return m.dydt[t, l] == m.dydl2[t, l] >>> model.pde = pyo.Constraint(model.t, model.l, rule=_pde_rule) Declare other model components and apply a discretization transformation ... Deactivate the differential equations at certain boundary points >>> for con in model.ode[:, model.t.first()]: ... con.deactivate() >>> for con in model.pde[0, :]: ... con.deactivate() >>> for con in model.pde[:, model.l.first()]: ... con.deactivate() >>> for con in model.pde[:, model.l.last()]: ... con.deactivate() Solve the model ... Note If you intend to use the pyomo.DAE Simulator on your model then you must use constraint deactivation instead of constraint skipping in the differential equation rule. Declaring Integrals Warning The Integral component is still under development and considered a prototype. It currently includes only basic functionality for simple integrals. We welcome feedback on the interface and functionality but we do not recommend using it on general models. Instead, integrals should be reformulated as differential equations. class pyomo.dae.Integral(*args, **kwds)[source] Represents an integral over a continuous domain The Integral component can be used to represent an integral taken over the entire domain of a ContinuousSet. Once every ContinuousSet in a model has been discretized, any integrals in the model will be converted to algebraic equations using the trapezoid rule. Future development will include more sophisticated numerical integration methods. Parameters: *args – Every indexing set needed to evaluate the integral expression wrt (ContinuousSet) – The continuous domain over which the integral is being taken rule (function) – Function returning the expression being integrated get_continuousset()[source] Return the ContinuousSet the integral is being taken over Declaring an Integral component is similar to declaring an Expression component. A simple example is shown below: >>> model = pyo.ConcreteModel() >>> model.time = ContinuousSet(bounds=(0,10)) >>> model.X = pyo.Var(model.time) >>> model.scale = pyo.Param(initialize=1E-3) >>> def _intX(m,t): ... return m.X[t] >>> model.intX = Integral(model.time,wrt=model.time,rule=_intX) >>> def _obj(m): ... return m.scale*m.intX >>> model.obj = pyo.Objective(rule=_obj) Notice that the positional arguments supplied to the Integral declaration must include all indices needed to evaluate the integral expression. The integral expression is defined in a function and supplied to the ‘rule’ keyword argument. Finally, a user must specify a ContinuousSet that the integral is being evaluated over. This is done using the ‘wrt’ keyword argument. Note The ContinuousSet specified using the ‘wrt’ keyword argument must be explicitly specified as one of the indexing sets (meaning it must be supplied as a positional argument). This is to ensure consistency in the ordering and dimension of the indexing sets After an Integral has been declared, it can be used just like a Pyomo Expression component and can be included in constraints or the objective function as shown above. If an Integral is specified with multiple positional arguments, i.e. multiple indexing sets, the final component will be indexed by all of those sets except for the ContinuousSet that the integral was taken over. In other words, the ContinuousSet specified with the ‘wrt’ keyword argument is removed from the indexing sets of the Integral even though it must be specified as a positional argument. This should become more clear with the following example showing a double integral over the ContinuousSet components model.t1 and model.t2. In addition, the expression is also indexed by the Set model.s. The mathematical representation and implementation in Pyomo are shown below: \[\sum_{s} \int_{t_2} \int_{t_1} \! X(t_1, t_2, s) \, dt_1 \, dt_2\] >>> model = pyo.ConcreteModel() >>> model.t1 = ContinuousSet(bounds=(0, 10)) >>> model.t2 = ContinuousSet(bounds=(-1, 1)) >>> model.s = pyo.Set(initialize=['A', 'B', 'C']) >>> model.X = pyo.Var(model.t1, model.t2, model.s) >>> def _intX1(m, t1, t2, s): ... return m.X[t1, t2, s] >>> model.intX1 = Integral(model.t1, model.t2, model.s, wrt=model.t1, ... rule=_intX1) >>> def _intX2(m, t2, s): ... return m.intX1[t2, s] >>> model.intX2 = Integral(model.t2, model.s, wrt=model.t2, rule=_intX2) >>> def _obj(m): ... return sum(m.intX2[k] for k in m.s) >>> model.obj = pyo.Objective(rule=_obj) Discretization Transformations Before a Pyomo model with DerivativeVar or Integral components can be sent to a solver it must first be sent through a discretization transformation. These transformations approximate any derivatives or integrals in the model by using a numerical method. The numerical methods currently included in pyomo.DAE discretize the continuous domains in the problem and introduce equality constraints which approximate the derivatives and integrals at the discretization points. Two families of discretization schemes have been implemented in pyomo.DAE, Finite Difference and Collocation. These schemes are described in more detail below. Note The schemes described here are for derivatives only. All integrals will be transformed using the trapezoid rule. The user must write a Python script in order to use these discretizations, they have not been tested on the pyomo command line. Example scripts are shown below for each of the discretization schemes. The transformations are applied to Pyomo model objects which can be further manipulated before being sent to a solver. Examples of this are also shown below. Finite Difference Transformation This transformation includes implementations of several finite difference methods. For example, the Backward Difference method (also called Implicit or Backward Euler) has been implemented. The discretization equations for this method are shown below: \[\begin{array}{l} \mathrm{Given: } \\ \frac{dx}{dt} = f(t, x) , \quad x(t_0) = x_{0} \\ \text{discretize $t$ and $x$ such that } \\ x(t_0 + kh) = x_{k} \\ x_{k + 1} = x_{k} + h * f(t_{k + 1}, x_{k + 1}) \\ t_{k + 1} = t_{k} + h \end{array}\]where \(h\) is the step size between discretization points or the size of each finite element. These equations are generated automatically as Constraints when the backward difference method is applied to a Pyomo model. There are several discretization options available to a dae.finite_difference transformation which can be specified as keyword arguments to the .apply_to() function of the transformation object. These keywords are summarized below: Keyword arguments for applying a finite difference transformation: ‘nfe’The desired number of finite element points to be included in the discretization. The default value is 10. ‘wrt’Indicates which ContinuousSet the transformation should be applied to. If this keyword argument is not specified then the same scheme will be applied to every ContinuousSet . ‘scheme’Indicates which finite difference method to apply. Options are ‘BACKWARD’, ‘CENTRAL’, or ‘FORWARD’. The default scheme is the backward difference method. If the existing number of finite element points in a ContinuousSet is less than the desired number, new discretization points will be added to the set. If a user specifies a number of finite element points which is less than the number of points already included in the ContinuousSet then the transformation will ignore the specified number and proceed with the larger set of points. Discretization points will never be removed from a ContinuousSet during the discretization. The following code is a Python script applying the backward difference method. The code also shows how to add a constraint to a discretized model. Discretize model using Backward Difference method >>> discretizer = pyo.TransformationFactory('dae.finite_difference') >>> discretizer.apply_to(model,nfe=20,wrt=model.time,scheme='BACKWARD') Add another constraint to discretized model >>> def _sum_limit(m): ... return sum(m.x1[i] for i in m.time) <= 50 >>> model.con_sum_limit = pyo.Constraint(rule=_sum_limit) Solve discretized model >>> solver = pyo.SolverFactory('ipopt') >>> results = solver.solve(model) Collocation Transformation This transformation uses orthogonal collocation to discretize the differential equations in the model. Currently, two types of collocation have been implemented. They both use Lagrange polynomials with either Gauss-Radau roots or Gauss-Legendre roots. For more information on orthogonal collocation and the discretization equations associated with this method please see chapter 10 of the book “Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes” by L.T. Biegler. The discretization options available to a dae.collocation transformation are the same as those described above for the finite difference transformation with different available schemes and the addition of the ‘ncp’ option. Additional keyword arguments for collocation discretizations: ‘scheme’The desired collocation scheme, either ‘LAGRANGE-RADAU’ or ‘LAGRANGE-LEGENDRE’. The default is ‘LAGRANGE-RADAU’. ‘ncp’The number of collocation points within each finite element. The default value is 3. Note If the user’s version of Python has access to the package Numpy then any number of collocation points may be specified, otherwise the maximum number is 10. Note Any points that exist in a ContinuousSet before discretization will be used as finite element boundaries and not as collocation points. The locations of the collocation points cannot be specified by the user, they must be generated by the transformation. The following code is a Python script applying collocation with Lagrange polynomials and Radau roots. The code also shows how to add an objective function to a discretized model. Discretize model using Radau Collocation >>> discretizer = pyo.TransformationFactory('dae.collocation') >>> discretizer.apply_to(model,nfe=20,ncp=6,scheme='LAGRANGE-RADAU') Add objective function after model has been discretized >>> def obj_rule(m): ... return sum((m.x[i]-m.x_ref)**2 for i in m.time) >>> model.obj = pyo.Objective(rule=obj_rule) Solve discretized model >>> solver = pyo.SolverFactory('ipopt') >>> results = solver.solve(model) Restricting Optimal Control Profiles When solving an optimal control problem a user may want to restrict the number of degrees of freedom for the control input by forcing, for example, a piecewise constant profile. Pyomo.DAE provides the reduce_collocation_points function to address this use-case. This function is used in conjunction with the dae.collocation discretization transformation to reduce the number of free collocation points within a finite element for a particular variable. class pyomo.dae.plugins.colloc.Collocation_Discretization_Transformation[source] reduce_collocation_points(instance, var=None, ncp=None, contset=None)[source] This method will add additional constraints to a model to reduce the number of free collocation points (degrees of freedom) for a particular variable. Parameters: instance (Pyomo model) – The discretized Pyomo model to add constraints to var (pyomo.environ.Var) – The Pyomo variable for which the degrees of freedom will be reduced ncp (int) – The new number of free collocation points for var. Must be less that the number of collocation points used in discretizing the model. contset (pyomo.dae.ContinuousSet) – The ContinuousSet that was discretized and for which the var will have a reduced number of degrees of freedom An example of using this function is shown below: >>> discretizer = pyo.TransformationFactory('dae.collocation') >>> discretizer.apply_to(model, nfe=10, ncp=6) >>> model = discretizer.reduce_collocation_points(model, ... var=model.u, ... ncp=1, ... contset=model.time) In the above example, the reduce_collocation_points function restricts the variable model.u to have only 1 free collocation point per finite element, thereby enforcing a piecewise constant profile. Fig. 1 shows the solution profile before and after applying the reduce_collocation_points function. Fig. 1 (left) Profile before applying the reduce_collocation_points function (right) Profile after applying the function, restricting model.u to have a piecewise constant profile. Applying Multiple Discretization Transformations Discretizations can be applied independently to each ContinuousSet in a model. This allows the user great flexibility in discretizing their model. For example the same numerical method can be applied with different resolutions: >>> discretizer = pyo.TransformationFactory('dae.finite_difference') >>> discretizer.apply_to(model,wrt=model.t1,nfe=10) >>> discretizer.apply_to(model,wrt=model.t2,nfe=100) This also allows the user to combine different methods. For example, applying the forward difference method to one ContinuousSet and the central finite difference method to another ContinuousSet: >>> discretizer = pyo.TransformationFactory('dae.finite_difference') >>> discretizer.apply_to(model,wrt=model.t1,scheme='FORWARD') >>> discretizer.apply_to(model,wrt=model.t2,scheme='CENTRAL') In addition, the user may combine finite difference and collocation discretizations. For example: >>> disc_fe = pyo.TransformationFactory('dae.finite_difference') >>> disc_fe.apply_to(model,wrt=model.t1,nfe=10) >>> disc_col = pyo.TransformationFactory('dae.collocation') >>> disc_col.apply_to(model,wrt=model.t2,nfe=10,ncp=5) If the user would like to apply the same discretization to all ContinuousSet components in a model, just specify the discretization once without the ‘wrt’ keyword argument. This will apply that scheme to all ContinuousSet components in the model that haven’t already been discretized. Custom Discretization Schemes A transformation framework along with certain utility functions has been created so that advanced users may easily implement custom discretization schemes other than those listed above. The transformation framework consists of the following steps: Specify Discretization Options Discretize the ContinuousSet(s) Update Model Components Add Discretization Equations Return Discretized Model If a user would like to create a custom finite difference scheme then they only have to worry about step (4) in the framework. The discretization equations for a particular scheme have been isolated from of the rest of the code for implementing the transformation. The function containing these discretization equations can be found at the top of the source code file for the transformation. For example, below is the function for the forward difference method: def _forward_transform(v,s): """ Applies the Forward Difference formula of order O(h) for first derivatives """ def _fwd_fun(i): tmp = sorted(s) idx = tmp.index(i) return 1/(tmp[idx+1]-tmp[idx])*(v(tmp[idx+1])-v(tmp[idx])) return _fwd_fun In this function, ‘v’ represents the continuous variable or function that the method is being applied to. ‘s’ represents the set of discrete points in the continuous domain. In order to implement a custom finite difference method, a user would have to copy the above function and just replace the equation next to the first return statement with their method. After implementing a custom finite difference method using the above function template, the only other change that must be made is to add the custom method to the ‘all_schemes’ dictionary in the dae.finite_difference class. In the case of a custom collocation method, changes will have to be made in steps (2) and (4) of the transformation framework. In addition to implementing the discretization equations, the user would also have to ensure that the desired collocation points are added to the ContinuousSet being discretized. Dynamic Model Simulation The pyomo.dae Simulator class can be used to simulate systems of ODEs and DAEs. It provides an interface to integrators available in other Python packages. Note The pyomo.dae Simulator does not include integrators directly. The user must have at least one of the supported Python packages installed in order to use this class. class pyomo.dae.Simulator(m, package='scipy')[source] Simulator objects allow a user to simulate a dynamic model formulated using pyomo.dae. Parameters: m (Pyomo Model) – The Pyomo model to be simulated should be passed as the first argument package (string) – The Python simulator package to use. Currently ‘scipy’ and ‘casadi’ are the only supported packages get_variable_order(vartype=None)[source] This function returns the ordered list of differential variable names. The order corresponds to the order being sent to the integrator function. Knowing the order allows users to provide initial conditions for the differential equations using a list or map the profiles returned by the simulate function to the Pyomo variables. Parameters: vartype (string or None) – Optional argument for specifying the type of variables to return the order for. The default behavior is to return the order of the differential variables. ‘time-varying’ will return the order of all the time-dependent algebraic variables identified in the model. ‘algebraic’ will return the order of algebraic variables used in the most recent call to the simulate function. ‘input’ will return the order of the time-dependent algebraic variables that were treated as inputs in the most recent call to the simulate function. Return type: list initialize_model()[source] This function will initialize the model using the profile obtained from simulating the dynamic model. simulate(numpoints=None, tstep=None, integrator=None, varying_inputs=None, initcon=None, integrator_options=None)[source] Simulate the model. Integrator-specific options may be specified as keyword arguments and will be passed on to the integrator. Parameters: numpoints (int) – The number of points for the profiles returned by the simulator. Default is 100 tstep (int or float) – The time step to use in the profiles returned by the simulator. This is not the time step used internally by the integrators. This is an optional parameter that may be specified in place of ‘numpoints’. integrator (string) – The string name of the integrator to use for simulation. The default is ‘lsoda’ when using Scipy and ‘idas’ when using CasADi varying_inputs (pyomo.environ.Suffix) – A Suffix object containing the piecewise constant profiles to be used for certain time-varying algebraic variables. initcon (list of floats) – The initial conditions for the the differential variables. This is an optional argument. If not specified then the simulator will use the current value of the differential variables at the lower bound of the ContinuousSet for the initial condition. integrator_options (dict) – Dictionary containing options that should be passed to the integrator. See the documentation for a specific integrator for a list of valid options. Returns: The first return value is a 1D array of time points corresponding to the second return value which is a 2D array of the profiles for the simulated differential and algebraic variables. Return type: numpy array, numpy array Note Any keyword options supported by the integrator may be specified as keyword options to the simulate function and will be passed to the integrator. Supported Simulator Packages The Simulator currently includes interfaces to SciPy and CasADi. ODE simulation is supported in both packages however, DAE simulation is only supported by CasADi. A list of available integrators for each package is given below. Please refer to the SciPy and CasADi documentation directly for the most up-to-date information about these packages and for more information about the various integrators and options. SciPy Integrators: ‘vode’ : Real-valued Variable-coefficient ODE solver, options for non-stiff and stiff systems ‘zvode’ : Complex-values Variable-coefficient ODE solver, options for non-stiff and stiff systems ‘lsoda’ : Real-values Variable-coefficient ODE solver, automatic switching of algorithms for non-stiff or stiff systems ‘dopri5’ : Explicit runge-kutta method of order (4)5 ODE solver ‘dop853’ : Explicit runge-kutta method of order 8(5,3) ODE solver CasADi Integrators: ‘cvodes’ : CVodes from the Sundials suite, solver for stiff or non-stiff ODE systems ‘idas’ : IDAS from the Sundials suite, DAE solver ‘collocation’ : Fixed-step implicit runge-kutta method, ODE/DAE solver ‘rk’ : Fixed-step explicit runge-kutta method, ODE solver Using the Simulator We now show how to use the Simulator to simulate the following system of ODEs: \[\begin{array}{l} \frac{d\theta}{dt} = \omega \\ \frac{d\omega}{dt} = -b*\omega -c*sin(\theta) \end{array}\]We begin by formulating the model using pyomo.DAE >>> m = pyo.ConcreteModel() >>> m.t = ContinuousSet(bounds=(0.0, 10.0)) >>> m.b = pyo.Param(initialize=0.25) >>> m.c = pyo.Param(initialize=5.0) >>> m.omega = pyo.Var(m.t) >>> m.theta = pyo.Var(m.t) >>> m.domegadt = DerivativeVar(m.omega, wrt=m.t) >>> m.dthetadt = DerivativeVar(m.theta, wrt=m.t) Setting the initial conditions >>> m.omega[0].fix(0.0) >>> m.theta[0].fix(3.14 - 0.1) >>> def _diffeq1(m, t): ... return m.domegadt[t] == -m.b * m.omega[t] - m.c * pyo.sin(m.theta[t]) >>> m.diffeq1 = pyo.Constraint(m.t, rule=_diffeq1) >>> def _diffeq2(m, t): ... return m.dthetadt[t] == m.omega[t] >>> m.diffeq2 = pyo.Constraint(m.t, rule=_diffeq2) Notice that the initial conditions are set by fixing the values of m.omega and m.theta at t=0 instead of being specified as extra equality constraints. Also notice that the differential equations are specified without using Constraint.Skip to skip enforcement at t=0. The Simulator cannot simulate any constraints that contain if-statements in their construction rules. To simulate the model you must first create a Simulator object. Building this object prepares the Pyomo model for simulation with a particular Python package and performs several checks on the model to ensure compatibility with the Simulator. Be sure to read through the list of limitations at the end of this section to understand the types of models supported by the Simulator. >>> sim = Simulator(m, package='scipy') After creating a Simulator object, the model can be simulated by calling the simulate function. Please see the API documentation for the Simulator for more information about the valid keyword arguments for this function. >>> tsim, profiles = sim.simulate(numpoints=100, integrator='vode') The simulate function returns numpy arrays containing time points and the corresponding values for the dynamic variable profiles. Simulator Limitations: Differential equations must be first-order and separable Model can only contain a single ContinuousSet Can’t simulate constraints with if-statements in the construction rules Need to provide initial conditions for dynamic states by setting the value or using fix() Specifying Time-Varying Inputs The Simulator supports simulation of a system of ODE’s or DAE’s with time-varying parameters or control inputs. Time-varying inputs can be specified using a Pyomo Suffix. We currently only support piecewise constant profiles. For more complex inputs defined by a continuous function of time we recommend adding an algebraic variable and constraint to your model. The profile for a time-varying input should be specified using a Python dictionary where the keys correspond to the switching times and the values correspond to the value of the input at a time point. A Suffix is then used to associate this dictionary with the appropriate Var or Param and pass the information to the Simulator. The code snippet below shows an example. >>> m = pyo.ConcreteModel() >>> m.t = ContinuousSet(bounds=(0.0, 20.0)) Time-varying inputs >>> m.b = pyo.Var(m.t) >>> m.c = pyo.Param(m.t, default=5.0) >>> m.omega = pyo.Var(m.t) >>> m.theta = pyo.Var(m.t) >>> m.domegadt = DerivativeVar(m.omega, wrt=m.t) >>> m.dthetadt = DerivativeVar(m.theta, wrt=m.t) Setting the initial conditions >>> m.omega[0] = 0.0 >>> m.theta[0] = 3.14 - 0.1 >>> def _diffeq1(m, t): ... return m.domegadt[t] == -m.b[t] * m.omega[t] - \ ... m.c[t] * pyo.sin(m.theta[t]) >>> m.diffeq1 = pyo.Constraint(m.t, rule=_diffeq1) >>> def _diffeq2(m, t): ... return m.dthetadt[t] == m.omega[t] >>> m.diffeq2 = pyo.Constraint(m.t, rule=_diffeq2) Specifying the piecewise constant inputs >>> b_profile = {0: 0.25, 15: 0.025} >>> c_profile = {0: 5.0, 7: 50} Declaring a Pyomo Suffix to pass the time-varying inputs to the Simulator >>> m.var_input = pyo.Suffix(direction=pyo.Suffix.LOCAL) >>> m.var_input[m.b] = b_profile >>> m.var_input[m.c] = c_profile Simulate the model using scipy >>> sim = Simulator(m, package='scipy') >>> tsim, profiles = sim.simulate(numpoints=100, ... integrator='vode', ... varying_inputs=m.var_input) Note The Simulator does not support multi-indexed inputs (i.e. if m.b in the above example was indexed by another set besides m.t) Dynamic Model Initialization Providing a good initial guess is an important factor in solving dynamic optimization problems. There are several model initialization tools under development in pyomo.DAE to help users initialize their models. These tools will be documented here as they become available. From Simulation The Simulator includes a function for initializing discretized dynamic optimization models using the profiles returned from the simulator. An example using this function is shown below Simulate the model using scipy >>> sim = Simulator(m, package='scipy') >>> tsim, profiles = sim.simulate(numpoints=100, integrator='vode', ... varying_inputs=m.var_input) Discretize the model using Orthogonal Collocation >>> discretizer = pyo.TransformationFactory('dae.collocation') >>> discretizer.apply_to(m, nfe=10, ncp=3) Initialize the discretized model using the simulator profiles >>> sim.initialize_model() Note A model must be simulated before it can be initialized using this function

```
pyomo.dae.ContinuousSet
```

**Pattern 6:** In addition, the user may combine finite difference and collocation discretizations. For example:

```
>>> disc_fe = pyo.TransformationFactory('dae.finite_difference')
>>> disc_fe.apply_to(model,wrt=model.t1,nfe=10)
>>> disc_col = pyo.TransformationFactory('dae.collocation')
>>> disc_col.apply_to(model,wrt=model.t2,nfe=10,ncp=5)
```

**Pattern 7:** Community Detection for Pyomo models This package separates model components (variables, constraints, and objectives) into different communities distinguished by the degree of connectivity between community members. Description of Package and detect_communities function The community detection package allows users to obtain a community map of a Pyomo model - a Python dictionary-like object that maps sequential integer values to communities within the Pyomo model. The package takes in a model, organizes the model components into a graph of nodes and edges, then uses Louvain community detection (Blondel et al, 2008) to determine the communities that exist within the model. In graph theory, a community is defined as a subset of nodes that have a greater degree of connectivity within themselves than they do with the rest of the nodes in the graph. In the context of Pyomo models, a community represents a subproblem within the overall optimization problem. Identifying these subproblems and then solving them independently can save computational work compared with trying to solve the entire model at once. Thus, it can be very useful to know the communities that exist in a model. The manner in which the graph of nodes and edges is constructed from the model directly affects the community detection. Thus, this package provides the user with a lot of control over the construction of the graph. The function we use for this community detection is shown below: pyomo.contrib.community_detection.detection.detect_communities(model, type_of_community_map='constraint', with_objective=True, weighted_graph=True, random_seed=None, use_only_active_components=True)[source] Detects communities in a Pyomo optimization model This function takes in a Pyomo optimization model and organizes the variables and constraints into a graph of nodes and edges. Then, by using Louvain community detection on the graph, a dictionary (community_map) is created, which maps (arbitrary) community keys to the detected communities within the model. Parameters: model (Block) – a Pyomo model or block to be used for community detection type_of_community_map (str, optional) – a string that specifies the type of community map to be returned, the default is ‘constraint’. ‘constraint’ returns a dictionary (community_map) with communities based on constraint nodes, ‘variable’ returns a dictionary (community_map) with communities based on variable nodes, ‘bipartite’ returns a dictionary (community_map) with communities based on a bipartite graph (both constraint and variable nodes) with_objective (bool, optional) – a Boolean argument that specifies whether or not the objective function is included in the model graph (and thus in ‘community_map’); the default is True weighted_graph (bool, optional) – a Boolean argument that specifies whether community_map is created based on a weighted model graph or an unweighted model graph; the default is True (type_of_community_map=’bipartite’ creates an unweighted model graph regardless of this parameter) random_seed (int, optional) – an integer that is used as the random seed for the (heuristic) Louvain community detection use_only_active_components (bool, optional) – a Boolean argument that specifies whether inactive constraints/objectives are included in the community map Returns: The CommunityMap object acts as a Python dictionary, mapping integer keys to tuples containing two lists (which contain the components in the given community) - a constraint list and variable list. Furthermore, the CommunityMap object stores relevant information about the given community map (dict), such as the model used to create it, its networkX representation, etc. Return type: CommunityMap object (dict-like object) As stated above, the characteristics of the NetworkX graph of the Pyomo model are very important to the community detection. The main graph features the user can specify are the type of community map, whether the graph is weighted or unweighted, and whether the objective function(s) is included in the graph generation. Below, the significance and reasoning behind including each of these options are explained in greater depth. Type of Community Map (type_of_community_map)In this package’s main function (detect_communities), the user can select 'bipartite', 'constraint', or 'variable' as an input for the ‘type_of_community_map’ argument, and these result in a community map based on a bipartite graph, a constraint node graph, or a variable node graph (respectively). If the user sets type_of_community_map='constraint', then each entry in the community map (which is a dictionary) contains a list of all the constraints in the community as well as all the variables contained in those constraints. For the model graph, a node is created for every active constraint in the model, an edge between two constraint nodes is created only if those two constraint equations share a variable, and the weight of each edge is equal to the number of variables the two constraint equations have in common. If the user sets type_of_community_map='variable', then each entry in the community map (which is a dictionary) contains a list of all the variables in the community as well as all the constraints that contain those variables. For the model graph, a node is created for every variable in the model, an edge between two variable nodes is created only if those two variables occur in the same constraint equation, and the weight of each edge is equal to the number of constraint equations in which the two variables occur together. If the user sets type_of_community_map='bipartite', then each entry in the community map (which is a dictionary) is simply all of the nodes in the community but split into a list of constraints and a list of variables. For the model graph, a node is created for every variable and every constraint in the model. An edge is created between a constraint node and a variable node only if the constraint equation contains the variable. (Edges are not drawn between nodes of the same type in a bipartite graph.) And as for the edge weights, the edges in the bipartite graph are unweighted regardless of what the user specifies for the weighted_graph parameter. (This is because for our purposes, the number of times a variable appears in a constraint is not particularly useful.) Weighted Graph/Unweighted Graph (weighted_graph)The Louvain community detection algorithm takes edge weights into account, so depending on whether the graph is weighted or unweighted, the communities that are found will vary. This can be valuable depending on how the user intends to use the community detection information. For example, if a user plans on feeding that information into an algorithm, the algorithm may be better suited to the communities detected in a weighted graph (or vice versa). With/Without Objective in the Graph (with_objective)This argument determines whether the objective function(s) will be included when creating the graphical representation of the model and thus whether the objective function(s) will be included in the community map. Some models have an objective function that contains so many of the model variables that it obscures potential communities within a model. Thus, it can be useful to call detect_communities(model, with_objective=False) on such a model to see whether isolating the other components of the model provides any new insights. External Packages NetworkX Python-Louvain The community detection package relies on two external packages, the NetworkX package and the Louvain community detection package. Both of these packages can be installed at the following URLs (respectively): https://pypi.org/project/networkx/ https://pypi.org/project/python-louvain/ The pip install and conda install commands are included below as well: pip install networkx pip install python-louvain conda install -c anaconda networkx conda install -c conda-forge python-louvain Usage Examples Let’s start off by taking a look at how we can use detect_communities to create a CommunityMap object. We’ll first use a model from Allman et al, 2019 : Required Imports >>> from pyomo.contrib.community_detection.detection import detect_communities, CommunityMap, generate_model_graph >>> from pyomo.contrib.mindtpy.tests.eight_process_problem import EightProcessFlowsheet >>> from pyomo.core import ConcreteModel, Var, Constraint >>> import networkx as nx Let's define a model for our use >>> def decode_model_1(): ... model = m = ConcreteModel() ... m.x1 = Var(initialize=-3) ... m.x2 = Var(initialize=-1) ... m.x3 = Var(initialize=-3) ... m.x4 = Var(initialize=-1) ... m.c1 = Constraint(expr=m.x1 + m.x2 <= 0) ... m.c2 = Constraint(expr=m.x1 - 3 * m.x2 <= 0) ... m.c3 = Constraint(expr=m.x2 + m.x3 + 4 * m.x4 ** 2 == 0) ... m.c4 = Constraint(expr=m.x3 + m.x4 <= 0) ... m.c5 = Constraint(expr=m.x3 ** 2 + m.x4 ** 2 - 10 == 0) ... return model >>> model = m = decode_model_1() >>> seed = 5 # To be used as a random seed value for the heuristic Louvain community detection Let's create an instance of the CommunityMap class (which is what gets returned by the function detect_communities): >>> community_map_object = detect_communities(model, type_of_community_map='bipartite', random_seed=seed) This community map object has many attributes that contain the relevant information about the community map itself (such as the parameters used to create it, the networkX representation, and other useful information). An important point to note is that the community_map attribute of the CommunityMap class is the actual dictionary that maps integers to the communities within the model. It is expected that the user will be most interested in the actual dictionary itself, so dict-like usage is permitted. If a user wishes to modify the actual dictionary (the community_map attribute of the CommunityMap object), creating a deep copy is highly recommended (or else any destructive modifications could have unintended consequences): new_community_map = copy.deepcopy(community_map_object.community_map) Let’s take a closer look at the actual community map object generated by detect_communities: >>> print(community_map_object) {0: (['c1', 'c2'], ['x1', 'x2']), 1: (['c3', 'c4', 'c5'], ['x3', 'x4'])} Printing a community map object is made to be user-friendly (by showing the community map with components replaced by their strings). However, if the default Pyomo representation of components is desired, then the community_map attribute or the repr() function can be used: >>> print(community_map_object.community_map) {0: ([<pyomo.core.base.constraint.ScalarConstraint object at ...>, <pyomo.core.base.constraint.ScalarConstraint object at ...>], [<pyomo.core.base.var.ScalarVar object at ...>, <pyomo.core.base.var.ScalarVar object at ...>]), 1: ([<pyomo.core.base.constraint.ScalarConstraint object at ...>, <pyomo.core.base.constraint.ScalarConstraint object at ...>, <pyomo.core.base.constraint.ScalarConstraint object at ...>], [<pyomo.core.base.var.ScalarVar object at ...>, <pyomo.core.base.var.ScalarVar object at ...>])} >>> print(repr(community_map_object)) {0: ([<pyomo.core.base.constraint.ScalarConstraint object at ...>, <pyomo.core.base.constraint.ScalarConstraint object at ...>], [<pyomo.core.base.var.ScalarVar object at ...>, <pyomo.core.base.var.ScalarVar object at ...>]), 1: ([<pyomo.core.base.constraint.ScalarConstraint object at ...>, <pyomo.core.base.constraint.ScalarConstraint object at ...>, <pyomo.core.base.constraint.ScalarConstraint object at ...>], [<pyomo.core.base.var.ScalarVar object at ...>, <pyomo.core.base.var.ScalarVar object at ...>])} generate_structured_model method of CommunityMap objectsIt may be useful to create a new model based on the communities found in the model - we can use the generate_structured_model method of the CommunityMap class to do this. Calling this method on a CommunityMap object returns a new model made up of blocks that correspond to each of the communities found in the original model. Let’s take a look at the example below: Use the CommunityMap object made from the first code example >>> structured_model = community_map_object.generate_structured_model() >>> structured_model.pprint() 2 Set Declarations b_index : Size=1, Index=None, Ordered=Insertion Key : Dimen : Domain : Size : Members None : 1 : Any : 2 : {0, 1} equality_constraint_list_index : Size=1, Index=None, Ordered=Insertion Key : Dimen : Domain : Size : Members None : 1 : Any : 1 : {1,} 1 Var Declarations x2 : Size=1, Index=None Key : Lower : Value : Upper : Fixed : Stale : Domain None : None : None : None : False : True : Reals 1 Constraint Declarations equality_constraint_list : Equality Constraints for the different forms of a given variable Size=1, Index=equality_constraint_list_index, Active=True Key : Lower : Body : Upper : Active 1 : 0.0 : b[0].x2 - x2 : 0.0 : True 1 Block Declarations b : Size=2, Index=b_index, Active=True b[0] : Active=True 2 Var Declarations x1 : Size=1, Index=None Key : Lower : Value : Upper : Fixed : Stale : Domain None : None : None : None : False : True : Reals x2 : Size=1, Index=None Key : Lower : Value : Upper : Fixed : Stale : Domain None : None : None : None : False : True : Reals 2 Constraint Declarations c1 : Size=1, Index=None, Active=True Key : Lower : Body : Upper : Active None : -Inf : b[0].x1 + b[0].x2 : 0.0 : True c2 : Size=1, Index=None, Active=True Key : Lower : Body : Upper : Active None : -Inf : b[0].x1 - 3*b[0].x2 : 0.0 : True 4 Declarations: x1 x2 c1 c2 b[1] : Active=True 2 Var Declarations x3 : Size=1, Index=None Key : Lower : Value : Upper : Fixed : Stale : Domain None : None : None : None : False : True : Reals x4 : Size=1, Index=None Key : Lower : Value : Upper : Fixed : Stale : Domain None : None : None : None : False : True : Reals 3 Constraint Declarations c3 : Size=1, Index=None, Active=True Key : Lower : Body : Upper : Active None : 0.0 : x2 + b[1].x3 + 4*b[1].x4**2 : 0.0 : True c4 : Size=1, Index=None, Active=True Key : Lower : Body : Upper : Active None : -Inf : b[1].x3 + b[1].x4 : 0.0 : True c5 : Size=1, Index=None, Active=True Key : Lower : Body : Upper : Active None : 0.0 : b[1].x3**2 + b[1].x4**2 - 10 : 0.0 : True 5 Declarations: x3 x4 c3 c4 c5 5 Declarations: b_index b x2 equality_constraint_list_index equality_constraint_list We see that there is an equality constraint list (equality_constraint_list) that has been created. This is due to the fact that the detect_communities function can return a community map that has Pyomo components (variables, constraints, or objectives) in more than one community, and thus, an equality_constraint_list is created to ensure that the new model still corresponds to the original model. This is explained in more detail below. Consider the case where community detection is done on a constraint node graph - this would result in communities that are made up of the corresponding constraints as well as all the variables that occur in the given constraints. Thus, it is possible for certain Pyomo components to be in multiple communities (and a similar argument exists for community detection done on a variable node graph). As a result, our structured model (the model returned by the generate_structured_model method) may need to have several “copies” of a certain component. For example, a variable original_model.x1 that exists in the original model may have corresponding forms structured_model.b[0].x1, structured_model.b[0].x1, structured_model.x1. In order for these components to meaningfully correspond to their counterparts in the original model, they must be bounded by equality constraints. Thus, we use an equality_constraint_list to bind different forms of a component from the original model. The last point to make about this method is that variables will be created outside of blocks if (1) an objective is not inside a block (for example if the community detection is done with_objective=False) or if (2) an objective/constraint contains a variable that is not in the same block as the given objective/constraint. visualize_model_graph method of CommunityMap objectsIf we want a visualization of the communities within the Pyomo model, we can use visualize_model_graph to do so. Let’s take a look at how this can be done in the following example: Create a CommunityMap object (so we can demonstrate the visualize_model_graph method) >>> community_map_object = cmo = detect_communities(model, type_of_community_map='bipartite', random_seed=seed) Generate a matplotlib figure (left_figure) - a constraint graph of the community map >>> left_figure, _ = cmo.visualize_model_graph(type_of_graph='constraint') >>> left_figure.show() Now, we will generate the figure on the right (a bipartite graph of the community map) >>> right_figure, _ = cmo.visualize_model_graph(type_of_graph='bipartite') >>> right_figure.show() An example of the two separate graphs created for these two function calls is shown below: These graph drawings very clearly demonstrate the communities within this model. The constraint graph (which is colored using the bipartite community map) shows a very simple illustration - one node for each constraint, with only one edge connecting the two communities (which represents the variable m.x2 common to m.c2 and m.c3 in separate communities) The bipartite graph is slightly more complicated and we can see again how there is only one edge between the two communities and more edges within each community. This is an ideal situation for breaking a model into separate communities since there is little connectivity between the communities. Also, note that we can choose different graph types (such as a variable node graph, constraint node graph, or bipartite graph) for a given community map. Let’s try a more complicated model (taken from Duran & Grossmann, 1986) - this example will demonstrate how the same graph can be illustrated using different community maps (in the previous example we illustrated different graphs with a single community map): Define the model >>> model = EightProcessFlowsheet() Now, we follow steps similar to the example above (see above for explanations) >>> community_map_object = cmo = detect_communities(model, type_of_community_map='constraint', random_seed=seed) >>> left_fig, pos = cmo.visualize_model_graph(type_of_graph='variable') >>> left_fig.show() As we did before, we will use the returned 'pos' to create a consistent graph layout >>> community_map_object = cmo = detect_communities(model, type_of_community_map='bipartite') >>> middle_fig, _ = cmo.visualize_model_graph(type_of_graph='variable', pos=pos) >>> middle_fig.show() >>> community_map_object = cmo = detect_communities(model, type_of_community_map='variable') >>> right_fig, _ = cmo.visualize_model_graph(type_of_graph='variable', pos=pos) >>> right_fig.show() We can see an example for the three separate graphs created by these three function calls below: The three graphs above are all variable graphs - which means the nodes represent variables in the model, and the edges represent constraint equations. The coloring differs because the three graphs rely on community maps that were created based on a constraint node graph, a bipartite graph, and a variable node graph (from left to right). For example, the community map that was generated from a constraint node graph (type_of_community_map='constraint') resulted in three communities (as seen by the purple, yellow, and blue nodes). generate_model_graph functionNow, we will take a look at generate_model_graph - this function can be used to create a NetworkX graph for a Pyomo model (and is used in detect_communities). Here, we will create a NetworkX graph from the model in our first example and then create the edge and adjacency list for the graph. generate_model_graph returns three things: a NetworkX graph of the given model a dictionary that maps the numbers used to represent the model components to the actual components (because Pyomo components cannot be directly added to a NetworkX graph) a dictionary that maps constraints to the variables in them. For this example, we will only need the NetworkX graph of the model and the number-to-component mapping. Define the model >>> model = decode_model_1() See above for the description of the items returned by 'generate_model_graph' >>> model_graph, number_component_map, constr_var_map = generate_model_graph(model, type_of_graph='constraint') The next two lines create and implement a mapping to change the node values from numbers into strings. The second line uses this mapping to create string_model_graph, which has the relabeled nodes (strings instead of numbers). >>> string_map = dict((number, str(comp)) for number, comp in number_component_map.items()) >>> string_model_graph = nx.relabel_nodes(model_graph, string_map) Now, we print the edge list and the adjacency list: Edge List: >>> for line in nx.generate_edgelist(string_model_graph): print(line) c1 c2 {'weight': 2} c1 c3 {'weight': 1} c2 c3 {'weight': 1} c3 c5 {'weight': 2} c3 c4 {'weight': 2} c4 c5 {'weight': 2} Adjacency List: >>> print(list(nx.generate_adjlist(string_model_graph))) ['c1 c2 c3', 'c2 c3', 'c3 c5 c4', 'c4 c5', 'c5'] It’s worth mentioning that in the code above, we do not have to create string_map to create an edge list or adjacency list, but for the sake of having an easily understandable output, it is quite helpful. (Without relabeling the nodes, the output below would not have the strings of the components but instead would have integer values.) This code will hopefully make it easier for a user to do the same. Functions in this Package Main module for community detection integration with Pyomo models. This module separates model components (variables, constraints, and objectives) into different communities distinguished by the degree of connectivity between community members. Original implementation developed by Rahul Joglekar in the Grossmann research group. class pyomo.contrib.community_detection.detection.CommunityMap(community_map, type_of_community_map, with_objective, weighted_graph, random_seed, use_only_active_components, model, graph, graph_node_mapping, constraint_variable_map, graph_partition)[source] This class is used to create CommunityMap objects which are returned by the detect_communities function. Instances of this class allow dict-like usage and store relevant information about the given community map, such as the model used to create them, their networkX representation, etc. The CommunityMap object acts as a Python dictionary, mapping integer keys to tuples containing two lists (which contain the components in the given community) - a constraint list and variable list. Methods: generate_structured_model visualize_model_graph generate_structured_model()[source] Using the community map and the original model used to create this community map, we will create structured_model, which will be based on the original model but will place variables, constraints, and objectives into or outside of various blocks (communities) based on the community map. Returns: structured_model – a Pyomo model that reflects the nature of the community map Return type: Block visualize_model_graph(type_of_graph='constraint', filename=None, pos=None)[source] This function draws a graph of the communities for a Pyomo model. The type_of_graph parameter is used to create either a variable-node graph, constraint-node graph, or bipartite graph of the Pyomo model. Then, the nodes are colored based on the communities they are in - which is based on the community map (self.community_map). A filename can be provided to save the figure, otherwise the figure is illustrated with matplotlib. Parameters: type_of_graph (str, optional) – a string that specifies the types of nodes drawn on the model graph, the default is ‘constraint’. ‘constraint’ draws a graph with constraint nodes, ‘variable’ draws a graph with variable nodes, ‘bipartite’ draws a bipartite graph (with both constraint and variable nodes) filename (str, optional) – a string that specifies a path for the model graph illustration to be saved pos (dict, optional) – a dictionary that maps node keys to their positions on the illustration Returns: fig (matplotlib figure) – the figure for the model graph drawing pos (dict) – a dictionary that maps node keys to their positions on the illustration - can be used to create consistent layouts for graphs of a given model pyomo.contrib.community_detection.detection.detect_communities(model, type_of_community_map='constraint', with_objective=True, weighted_graph=True, random_seed=None, use_only_active_components=True)[source] Detects communities in a Pyomo optimization model This function takes in a Pyomo optimization model and organizes the variables and constraints into a graph of nodes and edges. Then, by using Louvain community detection on the graph, a dictionary (community_map) is created, which maps (arbitrary) community keys to the detected communities within the model. Parameters: model (Block) – a Pyomo model or block to be used for community detection type_of_community_map (str, optional) – a string that specifies the type of community map to be returned, the default is ‘constraint’. ‘constraint’ returns a dictionary (community_map) with communities based on constraint nodes, ‘variable’ returns a dictionary (community_map) with communities based on variable nodes, ‘bipartite’ returns a dictionary (community_map) with communities based on a bipartite graph (both constraint and variable nodes) with_objective (bool, optional) – a Boolean argument that specifies whether or not the objective function is included in the model graph (and thus in ‘community_map’); the default is True weighted_graph (bool, optional) – a Boolean argument that specifies whether community_map is created based on a weighted model graph or an unweighted model graph; the default is True (type_of_community_map=’bipartite’ creates an unweighted model graph regardless of this parameter) random_seed (int, optional) – an integer that is used as the random seed for the (heuristic) Louvain community detection use_only_active_components (bool, optional) – a Boolean argument that specifies whether inactive constraints/objectives are included in the community map Returns: The CommunityMap object acts as a Python dictionary, mapping integer keys to tuples containing two lists (which contain the components in the given community) - a constraint list and variable list. Furthermore, the CommunityMap object stores relevant information about the given community map (dict), such as the model used to create it, its networkX representation, etc. Return type: CommunityMap object (dict-like object) Model Graph Generator Code pyomo.contrib.community_detection.community_graph.generate_model_graph(model, type_of_graph, with_objective=True, weighted_graph=True, use_only_active_components=True)[source] Creates a networkX graph of nodes and edges based on a Pyomo optimization model This function takes in a Pyomo optimization model, then creates a graphical representation of the model with specific features of the graph determined by the user (see Parameters below). (This function is designed to be called by detect_communities, but can be used solely for the purpose of creating model graphs as well.) Parameters: model (Block) – a Pyomo model or block to be used for community detection type_of_graph (str) – a string that specifies the type of graph that is created from the model ‘constraint’ creates a graph based on constraint nodes, ‘variable’ creates a graph based on variable nodes, ‘bipartite’ creates a graph based on constraint and variable nodes (bipartite graph). with_objective (bool, optional) – a Boolean argument that specifies whether or not the objective function is included in the graph; the default is True weighted_graph (bool, optional) – a Boolean argument that specifies whether a weighted or unweighted graph is to be created from the Pyomo model; the default is True (type_of_graph=’bipartite’ creates an unweighted graph regardless of this parameter) use_only_active_components (bool, optional) – a Boolean argument that specifies whether inactive constraints/objectives are included in the networkX graph Returns: bipartite_model_graph/projected_model_graph (nx.Graph) – a NetworkX graph with nodes and edges based on the given Pyomo optimization model number_component_map (dict) – a dictionary that (deterministically) maps a number to a component in the model constraint_variable_map (dict) – a dictionary that maps a numbered constraint to a list of (numbered) variables that appear in the constraint

```
detect_communities
```

**Pattern 8:** If we want a visualization of the communities within the Pyomo model, we can use visualize_model_graph to do so. Let’s take a look at how this can be done in the following example:

```
visualize_model_graph
```

## Reference Files

This skill includes comprehensive documentation in `references/`:

- **api.md** - Api documentation
- **explanation.md** - Explanation documentation
- **howto.md** - Howto documentation
- **index.html.md** - Index.Html documentation
- **modeling.md** - Modeling documentation
- **other.md** - Other documentation
- **reference.md** - Reference documentation
- **solvers.md** - Solvers documentation
- **tutorials.md** - Tutorials documentation

Use `view` to read specific reference files when detailed information is needed.

## Working with This Skill

### For Beginners
Start with the getting_started or tutorials reference files for foundational concepts.

### For Specific Features
Use the appropriate category reference file (api, guides, etc.) for detailed information.

### For Code Examples
The quick reference section above contains common patterns extracted from the official docs.

## Resources

### references/
Organized documentation extracted from official sources. These files contain:
- Detailed explanations
- Code examples with language annotations
- Links to original documentation
- Table of contents for quick navigation

### scripts/
Add helper scripts here for common automation tasks.

### assets/
Add templates, boilerplate, or example projects here.

## Notes

- This skill was automatically generated from official documentation
- Reference files preserve the structure and examples from source docs
- Code examples include language detection for better syntax highlighting
- Quick reference patterns are extracted from common usage examples in the docs

## Updating

To refresh this skill with updated documentation:
1. Re-run the scraper with the same configuration
2. The skill will be rebuilt with the latest information