Pyomo Network

Pyomo Network is a package that allows users to easily represent their model as a connected network of units. Units are blocks that contain ports, which contain variables, that are connected to other ports via arcs. The connection of two ports to each other via an arc typically represents a set of constraints equating each member of each port to each other, however there exist other connection rules as well, in addition to support for custom rules. Pyomo Network also includes a model transformation that will automatically expand the arcs and generate the appropriate constraints to produce an algebraic model that a solver can handle. Furthermore, the package also introduces a generic sequential decomposition tool that can leverage the modeling components to decompose a model and compute each unit in the model in a logically ordered sequence.

Modeling Components

Pyomo Network introduces two new modeling components to Pyomo:

pyomo.network.Port	A collection of variables, which may be connected to other ports
pyomo.network.Arc	Component used for connecting the members of two Port objects

Port

class pyomo.network.Port(*args, **kwds)

A collection of variables, which may be connected to other ports

The idea behind Ports is to create a bundle of variables that can be manipulated together by connecting them to other ports via Arcs. A preprocess transformation will look for Arcs and expand them into a series of constraints that involve the original variables contained within the Port. The way these constraints are built can be specified for each Port member when adding members to the port, but by default the Port members will be equated to each other. Additionally, other objects such as expressions can be added to Ports as long as they, or their indexed members, can be manipulated within constraint expressions.

Parameters:

• rule (function) – A function that returns a dict of (name: var) pairs to be initially added to the Port. Instead of var it could also be a tuples of (var, rule). Or it could return an iterable of either vars or tuples of (var, rule) for implied names.

• initialize – Follows same specifications as rule’s return value, gets initially added to the Port

• implicit – An iterable of names to be initially added to the Port as implicit vars

• extends (Port) – A Port whose vars will be added to this Port upon construction

static Equality(port, name, index_set)

Arc Expansion procedure to generate simple equality constraints

static Extensive(port, name, index_set, include_splitfrac=None, write_var_sum=True)

Arc Expansion procedure for extensive variable properties

This procedure is the rule to use when variable quantities should be conserved; that is, split for outlets and combined for inlets.

This will first go through every destination of the port (i.e., arcs whose source is this Port) and create a new variable on the arc’s expanded block of the same index as the current variable being processed to store the amount of the variable that flows over the arc. For ports that have multiple outgoing arcs, this procedure will create a single splitfrac variable on the arc’s expanded block as well. Then it will generate constraints for the new variable that relate it to the port member variable using the split fraction, ensuring that all extensive variables in the Port are split using the same ratio. The generation of the split fraction variable and constraint can be suppressed by setting the include_splitfrac argument to False.

Once all arc-specific variables are created, this procedure will create the “balancing constraint” that ensures that the sum of all the new variables equals the original port member variable. This constraint can be suppressed by setting the write_var_sum argument to False; in which case, a single constraint will be written that states the sum of the split fractions equals 1.

Finally, this procedure will go through every source for this port and create a new arc variable (unless it already exists), before generating the balancing constraint that ensures the sum of all the incoming new arc variables equals the original port variable.

Model simplifications:

If the port has a 1-to-1 connection on either side, it will not create the new variables and instead write a simple equality constraint for that side.

If the outlet side is not 1-to-1 but there is only one outlet, it will not create a splitfrac variable or write the split constraint, but it will still write the outsum constraint which will be a simple equality.

If the port only contains a single Extensive variable, the splitfrac variables and the splitting constraints will be skipped since they will be unnecessary. However, they can be still be included by passing include_splitfrac=True.

Note

If split fractions are skipped, the write_var_sum=False option is not allowed.

class pyomo.network.port._PortData(*args, **kwargs)

The following code snippet shows examples of declaring and using a Port component on a concrete Pyomo model:

>>> import pyomo.environ as pyo
>>> from pyomo.network import Port
>>> m = pyo.ConcreteModel()
>>> m.x = pyo.Var()
>>> m.y = pyo.Var(['a', 'b']) # can be indexed
>>> m.z = pyo.Var()
>>> m.e = 5 * m.z # you can add Pyomo expressions too
>>> m.w = pyo.Var()

>>> m.p = Port()
>>> m.p.add(m.x) # implicitly name the port member "x"
>>> m.p.add(m.y, "foo") # name the member "foo"
>>> m.p.add(m.e, rule=Port.Extensive) # specify a rule
>>> m.p.add(m.w, rule=Port.Extensive, write_var_sum=False) # keyword arg

Arc

class pyomo.network.Arc(*args, **kwds)

Component used for connecting the members of two Port objects

Parameters:

• source (Port) – A single Port for a directed arc. Aliases to src.

• destination (Port) – A single`Port for a directed arc. Aliases to dest.

• ports – A two-member list or tuple of single Ports for an undirected arc

• directed (bool) – Set True for directed. Use along with rule to be able to return an implied (source, destination) tuple.

• rule (function) – A function that returns either a dictionary of the arc arguments or a two-member iterable of ports

class pyomo.network.arc._ArcData(*args, **kwargs)

The following code snippet shows examples of declaring and using an Arc component on a concrete Pyomo model:

>>> import pyomo.environ as pyo
>>> from pyomo.network import Port, Arc
>>> m = pyo.ConcreteModel()
>>> m.x = pyo.Var()
>>> m.y = pyo.Var(['a', 'b'])
>>> m.u = pyo.Var()
>>> m.v = pyo.Var(['a', 'b'])
>>> m.w = pyo.Var()
>>> m.z = pyo.Var(['a', 'b']) # indexes need to match

>>> m.p = Port(initialize=[m.x, m.y])
>>> m.q = Port(initialize={"x": m.u, "y": m.v})
>>> m.r = Port(initialize={"x": m.w, "y": m.z}) # names need to match
>>> m.a = Arc(source=m.p, destination=m.q) # directed
>>> m.b = Arc(ports=(m.p, m.q)) # undirected
>>> m.c = Arc(ports=(m.p, m.q), directed=True) # directed
>>> m.d = Arc(src=m.p, dest=m.q) # aliases work
>>> m.e = Arc(source=m.r, dest=m.p) # ports can have both in and out

Arc Expansion Transformation

The examples above show how to declare and instantiate a Port and an Arc. These two components form the basis of the higher level representation of a connected network with sets of related variable quantities. Once a network model has been constructed, Pyomo Network implements a transformation that will expand all (active) arcs on the model and automatically generate the appropriate constraints. The constraints created for each port member will be indexed by the same indexing set as the port member itself.

During transformation, a new block is created on the model for each arc (located on the arc’s parent block), which serves to contain all of the auto generated constraints for that arc. At the end of the transformation, a reference is created on the arc that points to this new block, available via the arc property arc.expanded_block.

The constraints produced by this transformation depend on the rule assigned for each port member and can be different between members on the same port. For example, you can have two different members on a port where one member’s rule is Port.Equality and the other member’s rule is Port.Extensive.

Port.Equality is the default rule for port members. This rule simply generates equality constraints on the expanded block between the source port’s member and the destination port’s member. Another implemented expansion method is Port.Extensive, which essentially represents implied splitting and mixing of certain variable quantities. Users can refer to the documentation of the static method itself for more details on how this implicit splitting and mixing is implemented. Additionally, should users desire, the expansion API supports custom rules that can be implemented to generate whatever is needed for special cases.

The following code demonstrates how to call the transformation to expand the arcs on a model:

>>> import pyomo.environ as pyo
>>> from pyomo.network import Port, Arc
>>> m = pyo.ConcreteModel()
>>> m.x = pyo.Var()
>>> m.y = pyo.Var(['a', 'b'])
>>> m.u = pyo.Var()
>>> m.v = pyo.Var(['a', 'b'])

>>> m.p = Port(initialize=[m.x, (m.y, Port.Extensive)]) # rules must match
>>> m.q = Port(initialize={"x": m.u, "y": (m.v, Port.Extensive)})
>>> m.a = Arc(source=m.p, destination=m.q)

>>> pyo.TransformationFactory("network.expand_arcs").apply_to(m)

Sequential Decomposition

Pyomo Network implements a generic SequentialDecomposition tool that can be used to compute each unit in a network model in a logically ordered sequence.

The sequential decomposition procedure is commenced via the run method.

Creating a Graph

To begin this procedure, the Pyomo Network model is first utilized to create a networkx MultiDiGraph by adding edges to the graph for every arc on the model, where the nodes of the graph are the parent blocks of the source and destination ports. This is done via the create_graph method, which requires all arcs on the model to be both directed and already expanded. The MultiDiGraph class of networkx supports both direccted edges as well as having multiple edges between the same two nodes, so users can feel free to connect as many ports as desired between the same two units.

Computation Order

The order of computation is then determined by treating the resulting graph as a tree, starting at the roots of the tree, and making sure by the time each node is reached, all of its predecessors have already been computed. This is implemented through the calculation_order and tree_order methods. Before this, however, the procedure will first select a set of tear edges, if necessary, such that every loop in the graph is torn, while minimizing both the number of times any single loop is torn as well as the total number of tears.

Tear Selection

A set of tear edges can be selected in one of two ways. By default, a Pyomo MIP model is created and optimized resulting in an optimal set of tear edges. The implementation of this MIP model is based on a set of binary “torn” variables for every edge in the graph, and constraints on every loop in the graph that dictate that there must be at least one tear on the loop. Then there are two objectives (represented by a doubly weighted objective). The primary objective is to minimize the number of times any single loop is torn, and then secondary to that is to minimize the total number of tears. This process is implemented in the select_tear_mip method, which uses the model returned from the select_tear_mip_model method.

Alternatively, there is the select_tear_heuristic method. This uses a heuristic procedure that walks back and forth on the graph to find every optimal tear set, and returns each equally optimal tear set it finds. This method is much slower than the MIP method on larger models, but it maintains some use in the fact that it returns every possible optimal tear set.

A custom tear set can be assigned before calling the run method. This is useful so users can know what their tear set will be and thus what arcs will require guesses for uninitialized values. See the set_tear_set method for details.

Running the Sequential Decomposition Procedure

After all of this computational order preparation, the sequential decomposition procedure will then run through the graph in the order it has determined. Thus, the function that was passed to the run method will be called on every unit in sequence. This function can perform any arbitrary operations the user desires. The only thing that SequentialDecomposition expects from the function is that after returning from it, every variable on every outgoing port of the unit will be specified (i.e. it will have a set current value). Furthermore, the procedure guarantees to the user that for every unit, before the function is called, every variable on every incoming port of the unit will be fixed.

In between computing each of these units, port member values are passed across existing arcs involving the unit currently being computed. This means that after computing a unit, the expanded constraints from each arc coming out of this unit will be satisfied, and the values on the respective destination ports will be fixed at these new values. While running the computational order, values are not passed across tear edges, as tear edges represent locations in loops to stop computations (during iterations). This process continues until all units in the network have been computed. This concludes the “first pass run” of the network.

Guesses and Fixing Variables

When passing values across arcs while running the computational order, values at the destinations of each of these arcs will be fixed at the appropriate values. This is important to the fact that the procedure guarantees every inlet variable will be fixed before calling the function. However, since values are not passed across torn arcs, there is a need for user-supplied guesses for those values. See the set_guesses_for method for details on how to supply these values.

In addition to passing dictionaries of guesses for certain ports, users can also assign current values to the variables themselves and the procedure will pick these up and fix the variables in place. Alternatively, users can utilize the default_guess option to specify a value to use as a default guess for all free variables if they have no guess or current value. If a free variable has no guess or current value and there is no default guess option, then an error will be raised.

Similarly, if the procedure attempts to pass a value to a destination port member but that port member is already fixed and its fixed value is different from what is trying to be passed to it (by a tolerance specified by the almost_equal_tol option), then an error will be raised. Lastly, if there is more than one free variable in a constraint while trying to pass values across an arc, an error will be raised asking the user to fix more variables by the time values are passed across said arc.

Tear Convergence

After completing the first pass run of the network, the sequential decomposition procedure will proceed to converge all tear edges in the network (unless the user specifies not to, or if there are no tears). This process occurs separately for every strongly connected component (SCC) in the graph, and the SCCs are computed in a logical order such that each SCC is computed before other SCCs downstream of it (much like tree_order).

There are two implemented methods for converging tear edges: direct substitution and Wegstein acceleration. Both of these will iteratively run the computation order until every value in every tear arc has converged to within the specified tolerance. See the SequentialDecomposition parameter documentation for details on what can be controlled about this procedure.

The following code demonstrates basic usage of the SequentialDecomposition class:

>>> import pyomo.environ as pyo
>>> from pyomo.network import Port, Arc, SequentialDecomposition
>>> m = pyo.ConcreteModel()
>>> m.unit1 = pyo.Block()
>>> m.unit1.x = pyo.Var()
>>> m.unit1.y = pyo.Var(['a', 'b'])
>>> m.unit2 = pyo.Block()
>>> m.unit2.x = pyo.Var()
>>> m.unit2.y = pyo.Var(['a', 'b'])
>>> m.unit1.port = Port(initialize=[m.unit1.x, (m.unit1.y, Port.Extensive)])
>>> m.unit2.port = Port(initialize=[m.unit2.x, (m.unit2.y, Port.Extensive)])
>>> m.a = Arc(source=m.unit1.port, destination=m.unit2.port)
>>> pyo.TransformationFactory("network.expand_arcs").apply_to(m)

>>> m.unit1.x.fix(10)
>>> m.unit1.y['a'].fix(15)
>>> m.unit1.y['b'].fix(20)

>>> seq = SequentialDecomposition(tol=1.0E-3) # options can go to init
>>> seq.options.select_tear_method = "heuristic" # or set them like so
>>> # seq.set_tear_set([...]) # assign a custom tear set
>>> # seq.set_guesses_for(m.unit.inlet, {...}) # choose guesses
>>> def initialize(b):
...     # b.initialize()
...     pass
...
>>> seq.run(m, initialize)

class pyomo.network.SequentialDecomposition(**kwds)

A sequential decomposition tool for Pyomo Network models

The following parameters can be set upon construction of this class or via the options attribute.

Parameters:

• graph (MultiDiGraph) –

  A networkx graph representing the model to be solved.

  default=None (will compute it)

• tear_set (list) –

  A list of indexes representing edges to be torn. Can be set with a list of edge tuples via set_tear_set.

  default=None (will compute it)

• select_tear_method (str) –

  Which method to use to select a tear set, either “mip” or “heuristic”.

  default=”mip”

• run_first_pass (bool) –

  Boolean indicating whether or not to run through network before running the tear stream convergence procedure.

  default=True

• solve_tears (bool) –

  Boolean indicating whether or not to run iterations to converge tear streams.

  default=True

• guesses (ComponentMap) –

  ComponentMap of guesses to use for first pass (see set_guesses_for method).

  default=ComponentMap()

• default_guess (float) –

  Value to use if a free variable has no guess.

  default=None

• almost_equal_tol (float) –

  Difference below which numbers are considered equal when checking port value agreement.

  default=1.0E-8

• log_info (bool) –

  Set logger level to INFO during run.

  default=False

• tear_method (str) –

  Method to use for converging tear streams, either “Direct” or “Wegstein”.

  default=”Direct”

• iterLim (int) –

  Limit on the number of tear iterations.

  default=40

• tol (float) –

  Tolerance at which to stop tear iterations.

  default=1.0E-5

• tol_type (str) –

  Type of tolerance value, either “abs” (absolute) or “rel” (relative to current value).

  default=”abs”

• report_diffs (bool) –

  Report the matrix of differences across tear streams for every iteration.

  default=False

• accel_min (float) –

  Min value for Wegstein acceleration factor.

  default=-5

• accel_max (float) –

  Max value for Wegstein acceleration factor.

  default=0

• tear_solver (str) –

  Name of solver to use for select_tear_mip.

  default=”cplex”

• tear_solver_io (str) –

  Solver IO keyword for the above solver.

  default=None

• tear_solver_options (dict) –

  Keyword options to pass to solve method.

  default={}

calculation_order(G, roots=None, nodes=None)

Rely on tree_order to return a calculation order of nodes

Parameters:

• roots – List of nodes to consider as tree roots, if None then the actual roots are used

• nodes – Subset of nodes to consider in the tree, if None then all nodes are used

create_graph(model)

Returns a networkx MultiDiGraph of a Pyomo network model

The nodes are units and the edges follow Pyomo Arc objects. Nodes that get added to the graph are determined by the parent blocks of the source and destination Ports of every Arc in the model. Edges are added for each Arc using the direction specified by source and destination. All Arcs in the model will be used whether or not they are active (since this needs to be done after expansion), and they all need to be directed.

indexes_to_arcs(G, lst)

Converts a list of edge indexes to the corresponding Arcs

Parameters:

• G – A networkx graph corresponding to lst

• lst – A list of edge indexes to convert to tuples

Returns:

A list of arcs

run(model, function)

Compute a Pyomo Network model using sequential decomposition

Parameters:

• model – A Pyomo model

• function – A function to be called on each block/node in the network

select_tear_heuristic(G)

This finds optimal sets of tear edges based on two criteria. The primary objective is to minimize the maximum number of times any cycle is broken. The secondary criteria is to minimize the number of tears.

This function uses a branch and bound type approach.

Returns:

• tsets – List of lists of tear sets. All the tear sets returned are equally good. There are often a very large number of equally good tear sets.

• upperbound_loop – The max number of times any single loop is torn

• upperbound_total – The total number of loops

Improvements for the future

I think I can improve the efficiency of this, but it is good enough for now. Here are some ideas for improvement:

1. Reduce the number of redundant solutions. It is possible to find tears sets [1,2] and [2,1]. I eliminate redundant solutions from the results, but they can occur and it reduces efficiency.

2. Look at strongly connected components instead of whole graph. This would cut back on the size of graph we are looking at. The flowsheets are rarely one strongly connected component.

3. When you add an edge to a tear set you could reduce the size of the problem in the branch by only looking at strongly connected components with that edge removed.

4. This returns all equally good optimal tear sets. That may not really be necessary. For very large flowsheets, there could be an extremely large number of optimal tear edge sets.

select_tear_mip(G, solver, solver_io=None, solver_options={})

This finds optimal sets of tear edges based on two criteria. The primary objective is to minimize the maximum number of times any cycle is broken. The secondary criteria is to minimize the number of tears.

This function creates a MIP problem in Pyomo with a doubly weighted objective and solves it with the solver arguments.

select_tear_mip_model(G)

Generate a model for selecting tears from the given graph

Returns:

• model

• bin_list – A list of the binary variables representing each edge, indexed by the edge index of the graph

set_guesses_for(port, guesses)

Set the guesses for the given port

These guesses will be checked for all free variables that are encountered during the first pass run. If a free variable has no guess, its current value will be used. If its current value is None, the default_guess option will be used. If that is None, an error will be raised.

All port variables that are downstream of a non-tear edge will already be fixed. If there is a guess for a fixed variable, it will be silently ignored.

The guesses should be a dict that maps the following:

Port Member Name -> Value

Or, for indexed members, multiple dicts that map:

Port Member Name -> Index -> Value

For extensive members, “Value” must be a list of tuples of the form (arc, value) to guess a value for the expanded variable of the specified arc. However, if the arc connecting this port is a 1-to-1 arc with its peer, then there will be no expanded variable for the single arc, so a regular “Value” should be provided.

This dict cannot be used to pass guesses for variables within expression type members. Guesses for those variables must be assigned to the variable’s current value before calling run.

While this method makes things more convenient, all it does is:

self.options[“guesses”][port] = guesses

set_tear_set(tset)

Set a custom tear set to be used when running the decomposition

The procedure will use this custom tear set instead of finding its own, thus it can save some time. Additionally, this will be useful for knowing which edges will need guesses.

Parameters:

tset – A list of Arcs representing edges to tear

While this method makes things more convenient, all it does is:

self.options[“tear_set”] = tset

tear_set_arcs(G, method='mip', **kwds)

Call the specified tear selection method and return a list of arcs representing the selected tear edges.

The kwds will be passed to the method.

tree_order(adj, adjR, roots=None)

This function determines the ordering of nodes in a directed tree. This is a generic function that can operate on any given tree represented by the adjaceny and reverse adjacency lists. If the adjacency list does not represent a tree the results are not valid.

In the returned order, it is sometimes possible for more than one node to be calculated at once. So a list of lists is returned by this function. These represent a bredth first search order of the tree. Following the order, all nodes that lead to a particular node will be visited before it.

Parameters:

• adj – An adjeceny list for a directed tree. This uses generic integer node indexes, not node names from the graph itself. This allows this to be used on sub-graphs and graps of components more easily.

• adjR – The reverse adjacency list coresponing to adj

• roots – List of node indexes to start from. These do not need to be the root nodes of the tree, in some cases like when a node changes the changes may only affect nodes reachable in the tree from the changed node, in the case that roots are supplied not all the nodes in the tree may appear in the ordering. If no roots are supplied, the roots of the tree are used.