Structure-driven optimization methods for modular technical systems¶

The term structure-driven refers to the following idea. First, we partition the large, sparse nonlinear model of a technical system into smaller subproblems. Then, a suitable ordering of these smaller problems is determined so that the solution to the original problem can be reconstructed from the solutions of the smaller ones. The decomposition just sketched resembles the tree decomposition, which is very popular in discrete optimization and constraint satisfaction; this project deals with continuous problems. Evidence shows that appropriate decomposition can speed up the computations by several orders of magnitude.
Nonlinear programming is meant by optimization.
The image below shows an example of a modular technical system (a heating system, the image has been taken from the Modelica website). Component-based modeling is well suited for modular technical systems; Simulink and Dymola are examples of component-based modeling tools. The modules (components) of the technical systems help in decomposing the large model into smaller problems.

The source code is available on GitHub under the 3-clause BSD license.

THIS PROJECT IS A WORK IN PROGRESS.

Input¶

Currently, the models are written in AMPL, and after some black magic, the expressions are built-up in memory as a direct acyclic graph (DAG), using NetworkX DiGraph. It is somewhat similar to the expression trees in SymPy. Only nonlinear systems of equations are considered at the moment (steady-state modeling).

For example, the following AMPL code

var x; var y; var z;
equation: exp(3*x+2*y)+4*z = 1;

yields the directed acyclic graph below.

In the future, I would like to use either Python or JuMP for building the models. In any case, I will keep the AMPL interface too. Modelica is definitely on the agenda, accessed through the functional mock-up interface. This project is not aiming at creating yet another modeling environment: The goal is to plug the tools of this project into well-established modeling systems.

Reverse mode automatic differentiation¶

Source code is generated from the DAG representation of the expressions in order to compute the Jacobian with reverse mode automatic differentiation. Currently only Python code is emitted, in the near future, templated C++ code will also be generated. For example, for the above example exp(3*x+2*y)+4*z the following Python code is generated (hand-edited to improve readability):

# f = exp(3*x+2*y)+z
# Forward sweep
t1 = 3.0*x + 2.0*y
t2 = exp(t1)
f = 4.0*z + t2 - 1.0
# Backward sweep
u0 = 1.0
u1 = 4.0 * u0  # df/dz = 4
u2 = u0
u3 = t2 * u2
u4 = 3.0 * u3  # df/dx = 3*exp(3*x+2*y)
u5 = 2.0 * u3  # df/dy = 2*exp(3*x+2*y)

The templated C++ version of this code will greatly benefit from code optimization performed by the C++ compiler; I expect the generated code to be as good as hand-written.

Natural block structure¶

The modules of the technical systems partition the Jacobian into blocks in a fairly natural way. This natural block structure plays an important role in structure-driven algorithms: Computing the optimal partitioning and ordering of the smaller subproblems is NP-hard in the general case; one must resort to heuristics in practice. Independent results show that the heuristic which exploits the natural block structure often yields good quality partitioning and ordering.

The current way to pass the natural blocks is rather hackish: Suffixes are used, see Defining and using suffixes on page 302 in the AMPL book. In the future, component-based modeling tools will hopefully allow programmatic access to the natural block structure.

Minimum degree ordering¶

A basic minimum degree ordering has been implemented. The blue lines show the natural block structure.

My primary interest is chemical process modeling. The Jacobian of these models are very sparse but highly unsymmetric, numerically indefinite, not diagonally dominant and possibly ill-conditioned. There are many packages for the symmetric and slightly unsymmetric case. (In the slightly unsymmetric case, it is acceptable to introduce artificial fill-in to make the sparsity pattern symmetric and then use a sparse matrix ordering algorithm, developed for the symmetric case.) I have only found MC33 from the Harwell Subroutine Library that is applicable in the highly unsymmetric case. Since MC33 is based on a heuristic, it unfortunately fails on those chemical process models that are of interest to me.

Graph coloring¶

Depending on the implementation, efficient forward-mode automatic differentiation may require well-chosen seed vectors; these seeds can be computed with graph coloring. Even though graph coloring is NP-complete in general, the minimum degree ordering enables an efficient greedy coloring heuristic.