Imperial College London


Faculty of EngineeringDepartment of Chemical Engineering

Reader in Process Systems Engineering



b.chachuat Website




609Roderic Hill BuildingSouth Kensington Campus





My research interests are in the general area of Process Systems Engineering (PSE), with emphasis on global optimization, optimization-based process control, and numerical computation. My primary research objective is to develop new methods, concepts and tools for the fast and reliable optimization of complex process systems.

Applications in the field of process design and operations aim to develop safe, environmentally and economically sustainable (bio)chemical processes, including micro power generation, integrated biotechnological processes for bioenergy production and waste treatment. Other applications are in the field of signal transduction in systems biology.

Representative on-going projects in my group include:

Global Optimization of Transient Processes

Dynamic optimization problems encountered in various engineering and scientific fields frequently exhibit multiple suboptimal solutions. A prototypical example is in the fields of model identification, where failure to determine the best possible fit can lead to false conclusions regarding the validity of a candidate model. Other applications are in the field of robust and scenario-integration optimization.
Here, we are investigating deterministic global optimization methods, such as spatial branch-and-bound, which have the appealing property that a global optimum can be located, at an arbitrary precision, in a finite number of itterations. Pivotal to the success of these procedures is the ability to compute tight bounds or, better, convex/concave relaxations for the solutions of the underlying differential equations. Our current focus is on discretize-then-relax approaches that build upon verified ODE methods to compute the bounds, and their extension to handle more general DAE systems.

Dynamic Real-time Optimization (RTO) of Uncertain Processes

Challenges in real-time process optimization arise from the inability to build and adapt accurate models for complex systems, both steady-state and transient systems. In this research, we investigate different ways of using input/output data to compensate for model uncertainty and process disturbances, thereby promoting optimal operation. Our recent efforts have been directed towards the development of the modifier-adaptation and NCO-tracking technologies, in particular systematic design procedures. The combination of various (D)RTO paradigms along with their better integration/coordination with lower-level control systems are under investigation.

Design and Operation of Synthetic Ecosystems for Bioenergy Production and Waste Treatment

Microbial ecosystems, such as anaerobic digestion and activated sludge, are ubiquitous in today's society for the treatment of waste effluents. There exist strong environmental and economical incentives for designing and operating these biotechnological processes in the most sustainable way. The combination of several microbial ecosystems in an integrated manner, also known as synthetic ecosystems, holds much promise in this context. They open the possibility for treating multiple pollution sources and producing various forms of bioenergy, all in the same process.

The study of combined microalgal/bacterial ecosystems, through the coupling of a photobioreactor and an anaerobic digester, is currently underway in collaboration with LBE-INRA (Narbonne, France), and  COMORE-INRIA (Sophia-Antipolis, France). The foremost advantage of this integrated solution is the recovery nutrients—mainly inorganic nitrogen and phosphorus, along with some trace elements—from the biomass residue left behind after oil extraction from microalgae. Besides limiting the need for exogenous nutrient sources, raw or pretreated wastewater can also be fed to the anaerobic digester as a means to control the carbon-to-nitrogen ratio, while treating the organic pollution in the wastewater. From a PSE standpoint, such biotechnological processes exhibit significantly higher complexity, and the difficulty is exacerbated by the large amount of modeling uncertainty. On top of this, microbial ecosystems are inherently time-varying due to changing environmental conditions and microbial adaptation, and exhibit strong nonlinear phenomena that span multiple time scales, ranging from seconds to weeks. In this research, we aim to develop the methods and tools needed to bring this new generation of sustainable bioprocess to a reality and pave the road towards their successful deployment.


A Versatile Library for McCormick Relaxations and Taylor Models

MC++ is a library that computes convex/concave bounds and Taylor models for enclosing the range of factorable functions. Factorable functions are those defined as a finite recursive composition of binary sums, binary products and univariate functions; this is an extremely inclusive class of functions containing nearly every function which can be represented finitely on a computer. The main applications of MC++ are in the area of deterministic global optimization as well as for the verified solution of nonlinear algebraic equations and ordinary differential equations (ODEs). MC++ is programmed in C++ and makes extensive use of class templates and operator overloading. Although less performant than source code tranformation, this approach offers great flexibility and is particularly well suited for proof-of-concept implementations. Moreover, it is perfectly adequate for many small- to medium-size problems.

MC++ is the successor of libMC. It features two data types: (i) the class McCormick computes convex/concave bounds and subgradients for factorable functions, based on the generalized McCormick technique; and (ii) the class Taylor computes Taylor models for factorable functions as a means to bound their range, based on Taylor model arithmetic. Both McCormick and Taylor are template classes, which can be combined, for example by using McCormick inside Taylor and vice versa. Moreover, the template parameter can be used to pass interval data types, such as PROFIL, FILIB++ or INTLIB.

The data types McCormick and Taylor can also be used as template parameters in the library FADBAD++, that implements the forward, backward and Taylor methods of automatic differentiation (AD). This way, not only can bounds and relaxations be computed for factorable functions, but also bounds and relaxations for the derivatives and/or the Taylor coefficients of these functions.

McCormick Relaxations and Subgradient Propagation for Factorable Functions

Convex and concave relaxations provide an alternative, usually tighter, way of computing bounds on the range of functions than classical interval analysis. Given a convex set P in IRn and a function ƒ :P → IRn, a convex function ƒcv:P → IRn and a concave function ƒcc:P → IRn are called a convex relaxation and a concave relaxation of f on P, respectively, if ƒcv(p) ≤ ƒ(p) ≤ ƒcc(p) for all p in P.

In the adjacent figure, a function ƒ, a convex relaxation ƒcv and a concave relaxation ƒcc are represented in red, green and blue, respectively. Unlike the intervals bounds ƒL and ƒU, convex and concave relaxations are not constant on P. The values of convex and concave relaxations of ƒ at a given p in P are called convex/concave bounds at p and are denoted by [ƒcvcc](p).

For factorable functions, the construction of convex and concave bounds is possible through several mechanisms; MC++ implements the generalized McCormick technique and its extension for subgradient propagation.

• McCormick, G. P., Computability of Global Solutions to Factorable Nonconvex Programs: Part I–Convex Underestimating Problems, Mathematical Programming 10:147-175, 1976
• Mitsos, A., Chachuat, B., and Barton, P. I., McCormick-Based Relaxations of Algorithms, SIAM Journal on Optimization, 20(2):573-601, 2009
• Scott, J. K., Stuber, M. D., and Barton, P. I., Generalized McCormick Relaxations, Journal of Global Optimization, in press

Crash course in Computing Convex/concave Bounds and Subgradients

Suppose we are interested in computing convex/concave bounds for the real-valued function ƒ(x,y)=x exp(x+y2) - y2 for (x,y) in [1,2]x[0,1] at the point (x,y)=(1.5,0.5).

First, we need to specify the underlying data type for interval computations. The INTERVAL type in the PROFIL library is considered here:

#include "mccormick.h" #include "mcprofil.h"       typedef INTERVAL I;
      typedef mc::McCormick<I> MC;

Next, the variables x and y are defined as follows:

MC X( I(1., 2.), 1.5 ); MC Y( I(0., 1.), 0.5 );

Essentially, the first line means that X is a variable of class McCormick, that it belongs to the interval [1,2], and that its current value is 1.5. Likewise,  y is a variable of class McCormick that it belongs to the interval [1,2] and has current value 0.5.

Once x and y have been defined, the convex/concave bounds of ƒ on [1,2]x[0,1] at (1.5,0.5) are simply calculated as follows:

MC F = X * exp( X + pow( Y, 2 ) ) - pow( Y, 2 );

The values of the convex and concave bounds are retrieved as:

double Fcvx =;
      double Fccv =;

Likewise, lower and upper bounds on the range of ƒ on [1,2]x[0,1] are retreived as:

double Flb = F.l();
      double Fub = F.u();

The computation of a subgradient of a McCormick relaxation requires the specification of the independent variables via the .sub method, prior to the function evaluation. In the previous example, defining x and y as the subgradient components 0 and 1, respectively, is done as (note that the index starts at 0 not 1!):

X.sub( 2, 0 ); Y.sub( 2, 1 );

The convex and concave bounds and subgradients are calculated as previously, and the results can be retreived as follows:

MC F = X * exp( X + pow( Y, 2 ) ) - pow( Y, 2 );
      const double* dFcvx = F.cvsub();
      const double* dFccv = F.ccsub();

The convex/concave bounds and subgradients can also be displayed as:

std::cout << F;

By repeating the computation of convex/concave bounds at a large number of points in [1,2]x[0,1], the actual convex and concave relaxations can be generated and displayed. The McCormick relaxations of ƒ are shown in the figure below.

Taylor Model Computation for Factorable Functions

Taylor models provide yet another way of computing tight bounds on the range of functions. Initially, Taylor models were introduced to reduce the overestimation in interval analysis, by combining interval arithmetic with symbolic computations. Given a set P in IRn and a function ƒ :P → IRn, a qth-order Taylor model of ƒ on P is defined as Τƒ = Pƒ + Rƒ, where Pƒ is an np-variate polynomial of order q and Rƒ = [rƒL,rƒU] is an interval vector, such that ƒ(p) ε Pƒ(p) + Rƒ, for all p ε P. The range of a Taylor model Τƒ on P is the set {Pƒ(p) + Rƒ: p ε P}, which turns out to be an interval and requires bounding the range of the polynomial Pƒ.

In the adjacent figure, a univariate function ƒ and Taylor model Τƒ of that function on the interval P=[pL,pU] are depicted in red and green/blue colors, respectively. Observe, in particular, that Τƒ encloses ƒ between two hypersurfaces.
Similar to convex/concave bounds, Taylor models can be computed for any factorable function, provided that the participating univariate functions are q+1 times continuously differentiable on their domain.  The polynomial part is propagated by symbolic calculations wherever possible, the interval remainder term and all polynomial terms of order higher than q are processed according to the rules of interval arithmetic.

The implementation of Taylor models in MC++ is not validated in the sense that the round-off errors due to symbolic operations in floating-point arithmetic are not accounted for. Regarding range bounding, MC++ considers the first-order and diagonal second-order terms for exact bounding only, while the remaining terms are directly evaluated via interval analysis.

• Makino, K., and Berz, M., "Taylor models and other validated functional methods", Int J Pure Appl Math 4:379-456, 2003
• Neumaier, A., "Taylor forms - Use and limits", Reliable Computing, 9(1):43-79, 2002

Crash course in Computing Taylor Models

Suppose we are interested in computing a Taylor model for the real-valued function ƒ(x,y)=x exp(x+y2) - y2 for (x,y) in [1,2]x[0,1] at the point (x,y)=(1.5,0.5).

The INTERVAL type in the PROFIL library is used here again to perform the underlying interval computations:

#include "taylor.h" #include "mcprofil.h"       typedef INTERVAL I;
      typedef mc::Taylor<I> TM;

First, the number of independent variables (2 in this example) in the Taylor model as well as the order of the Taylor model (4th order here) are defined:

TM::size( 2, 4 );

Then, the variables x and y are defined as follows:

TM X( 0, I(1.,2.) ); TM Y( 1, I(0.,1.) );

Essentially, the first line means that X is a variable of class Taylor, that it belongs to the interval [1,2], and that it has index 0. The same holds for the McCormick variable Y that belongs to the interval [0,1] and has index 1.

Once x and y have been defined, a Taylor modelΤƒ of ƒ on [1,2]x[0,1] is simply calculated as follows:

TM F =  x * exp( X + pow( Y, 2 ) ) - pow( Y, 2 );

Taylor model derived bounds for the range of ƒ on [1,2]x[0,1] are then computed as:

I B = F.bound();

Detailed information on the Taylor model and the corresponding bounds can be displayed as:

std::cout << F;

The resulting 4th-order Taylor model of ƒ and the corresponding interval bounds are shown in the figure below.