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:
A collection of variables, which may be connected to other ports 

Component used for connecting the members of two Port objects 
Port
 class pyomo.network.Port(*args, **kwds)[source]
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)[source]
Arc Expansion procedure to generate simple equality constraints
 static Extensive(port, name, index_set, include_splitfrac=None, write_var_sum=True)[source]
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 arcspecific 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 1to1 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 1to1 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.
The following code snippet shows examples of declaring and using a
Port
component on a
concrete Pyomo model:
>>> from pyomo.environ import *
>>> from pyomo.network import *
>>> m = ConcreteModel()
>>> m.x = Var()
>>> m.y = Var(['a', 'b']) # can be indexed
>>> m.z = Var()
>>> m.e = 5 * m.z # you can add Pyomo expressions too
>>> m.w = 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)[source]
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 twomember 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 twomember iterable of ports
The following code snippet shows examples of declaring and using an
Arc
component on a
concrete Pyomo model:
>>> from pyomo.environ import *
>>> from pyomo.network import *
>>> m = ConcreteModel()
>>> m.x = Var()
>>> m.y = Var(['a', 'b'])
>>> m.u = Var()
>>> m.v = Var(['a', 'b'])
>>> m.w = Var()
>>> m.z = 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:
>>> from pyomo.environ import *
>>> from pyomo.network import *
>>> m = ConcreteModel()
>>> m.x = Var()
>>> m.y = Var(['a', 'b'])
>>> m.u = Var()
>>> m.v = 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)
>>> 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
usersupplied 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:
>>> from pyomo.environ import *
>>> from pyomo.network import *
>>> m = ConcreteModel()
>>> m.unit1 = Block()
>>> m.unit1.x = Var()
>>> m.unit1.y = Var(['a', 'b'])
>>> m.unit2 = Block()
>>> m.unit2.x = Var()
>>> m.unit2.y = 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)
>>> 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.0E3) # 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)[source]
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.0E8
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.0E5
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)[source]
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)[source]
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)[source]
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={})[source]
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)[source]
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)[source]
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 nontear 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 1to1 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)[source]
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)[source]
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 subgraphs 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.