Trust Region Framework Method Solver

The Trust Region Framework (TRF) method solver allows users to solve hybrid glass box/black box optimization problems in which parts of the system are modeled with open, equation-based models and parts of the system are black boxes. This method utilizes surrogate models that substitute high-fidelity models with low-fidelity basis functions, thus avoiding the direct implementation of the large, computationally expensive high-fidelity models. This is done iteratively, resulting in fewer calls to the computationally expensive functions.

This module implements the method from Yoshio & Biegler [Yoshio & Biegler, 2021] and represents a rewrite of the original 2018 implementation of the algorithm from Eason & Biegler [Eason & Biegler, 2018].

In the context of this updated module, black box functions are implemented as Pyomo External Functions.

This work was conducted as part of the Institute for the Design of Advanced Energy Systems (IDAES) with support through the Simulation-Based Engineering, Crosscutting Research Program within the U.S. Department of Energy’s Office of Fossil Energy and Carbon Management.

Methodology Overview

The formulation of the original hybrid problem is:

\[\begin{array}{ll} \min_{} & f\left(z, w, d\left(w\right)\right) \\ \text{s.t.} & h\left(z, w, d\left(w\right)\right) = 0 \\ & g\left(z, w, d\left(w\right)\right) \leq 0 \end{array}\]

where:

  • \(w \in \mathbb{R}^m\) are the inputs to the external functions

  • \(z \in \mathbb{R}^n\) are the remaining decision variables (i.e., degrees of freedom)

  • \(d(w) : \mathbb{R}^m \to \mathbb{R}^p\) are the outputs of the external functions as a function of \(w\)

  • \(f\), h, g, d are all assumed to be twice continuously differentiable

This formulation is reworked to separate all external function information as follows to enable the usage of the trust region method:

\[\begin{array}{ll} \min_{x} & f\left(x\right) \\ \text{s.t.} & h\left(x\right) = 0 \\ & g\left(x\right) \leq 0 \\ & y = d\left(w\right) \end{array}\]

where:

  • \(y \in \mathbb{R}^p\) are the outputs of the external functions

  • \(x^T = [w^T, y^T, z^T]\) is a set of all inputs and outputs

Using this formulation and a user-supplied low-fidelity/ideal model basis function \(b\left(w\right)\), the algorithm iteratively solves subproblems using the surrogate model:

\[r_k\left(w\right) = b\left(w\right) + \left( d\left(w_k\right) - b\left(w_k\right) \right) + \left( \nabla d\left(w_k\right) - \nabla b\left(w_k\right) \right)^T \left( w - w_k \right)\]

This acts similarly to Newton’s method in that small, incremental steps are taken towards an optimal solution. At each iteration, the current solution of the subproblem is compared to the previous solution to ensure that the iteration has moved in a direction towards an optimal solution. If not true, the step is rejected. If true, the step is accepted and the surrogate model is updated for the next iteration.

When using TRF, please consider citing the above papers.

TRF Inputs

The required inputs to the TRF solve method are the following:

  • The optimization model

  • List of degree of freedom variables within the model

The optional input to the TRF solve method is the following:

  • The external function surrogate model rule (“basis function”)

TRF Solver Interface

Note

The keyword arguments can be updated at solver instantiation or later when the solve method is called.

class pyomo.contrib.trustregion.TRF.TrustRegionSolver(**kwds)[source]

The Trust Region Solver is a ‘solver’ based on the 2016/2018/2020 AiChE papers by Eason (2016/2018), Yoshio (2020), and Biegler.

solve(model, degrees_of_freedom_variables, ext_fcn_surrogate_map_rule=None, **kwds)[source]

This method calls the TRF algorithm.

Parameters:

  • model (ConcreteModel) – The model to be solved using the Trust Region Framework.

  • degrees_of_freedom_variables (List[Var]) – User-supplied input. The user must provide a list of vars which are the degrees of freedom or decision variables within the model.

  • ext_fcn_surrogate_map_rule (Function, optional) – In the 2020 Yoshio/Biegler paper, this is referred to as the basis function b(w). This is the low-fidelity model with which to solve the original process model problem and which is integrated into the surrogate model. The default is 0 (i.e., no basis function rule.)

Keyword Arguments:

  • solver (default='ipopt') – Solver to use. Default = ipopt.

  • keepfiles (Bool, default=False) – Optional. Whether or not to write files of sub-problems for use in debugging. Default = False.

  • tee (Bool, default=False) – Optional. Sets the tee for sub-solver(s) utilized. Default = False.

  • verbose (Bool, default=False) – Optional. When True, print each iteration’s relevant information to the console as well as to the log. Default = False.

  • trust_radius (PositiveFloat, default=1.0) – Initial trust region radius delta_0. Default = 1.0.

  • minimum_radius (PositiveFloat, default=1e-06) – Minimum allowed trust region radius delta_min. Default = 1e-6.

  • maximum_radius (PositiveFloat, default=100.0) – Maximum allowed trust region radius. If trust region radius reaches maximum allowed, solver will exit. Default = 100 * trust_radius.

  • maximum_iterations (PositiveInt, default=50) – Maximum allowed number of iterations. Default = 50.

  • feasibility_termination (PositiveFloat, default=1e-05) – Feasibility measure termination tolerance epsilon_theta. Default = 1e-5.

  • step_size_termination (PositiveFloat, default=1e-05) – Step size termination tolerance epsilon_s. Matches the feasibility termination tolerance by default.

  • minimum_feasibility (PositiveFloat, default=0.0001) – Minimum feasibility measure theta_min. Default = 1e-4.

  • switch_condition_kappa_theta (In(0..1), default=0.1) – Switching condition parameter kappa_theta. Contained in open set (0, 1). Default = 0.1.

  • switch_condition_gamma_s (PositiveFloat, default=2.0) – Switching condition parameter gamma_s. Must satisfy: gamma_s > 1/(1+mu) where mu is contained in set (0, 1]. Default = 2.0.

  • radius_update_param_gamma_c (In(0..1), default=0.5) – Lower trust region update parameter gamma_c. Default = 0.5.

  • radius_update_param_gamma_e (In[1..inf], default=2.5) – Upper trust region update parameter gamma_e. Default = 2.5.

  • ratio_test_param_eta_1 (In(0..1), default=0.05) – Lower ratio test parameter eta_1. Must satisfy: 0 < eta_1 <= eta_2 < 1. Default = 0.05.

  • ratio_test_param_eta_2 (In(0..1), default=0.2) – Lower ratio test parameter eta_2. Must satisfy: 0 < eta_1 <= eta_2 < 1. Default = 0.2.

  • maximum_feasibility (PositiveFloat, default=50.0) – Maximum allowable feasibility measure theta_max. Parameter for use in filter method.Default = 50.0.

  • param_filter_gamma_theta (In(0..1), default=0.01) – Fixed filter parameter gamma_theta within (0, 1). Default = 0.01

  • param_filter_gamma_f (In(0..1), default=0.01) – Fixed filter parameter gamma_f within (0, 1). Default = 0.01

TRF Usage Example

Two examples can be found in the examples subdirectory. One of them is implemented below.

Step 0: Import Pyomo

>>> # === Required imports ===
>>> import pyomo.environ as pyo

Step 1: Define the external function and its gradient

>>> # === Define a 'black box' function and its gradient ===
>>> def ext_fcn(a, b):
...     return pyo.sin(a - b)
>>> def grad_ext_fcn(args, fixed):
...     a, b = args[:2]
...     return [ pyo.cos(a - b), -pyo.cos(a - b) ]

Step 2: Create the model

>>> # === Construct the Pyomo model object ===
>>> def create_model():
...     m = pyo.ConcreteModel()
...     m.name = 'Example 1: Eason'
...     m.z = pyo.Var(range(3), domain=pyo.Reals, initialize=2.)
...     m.x = pyo.Var(range(2), initialize=2.)
...     m.x[1] = 1.0
...
...     m.ext_fcn = pyo.ExternalFunction(ext_fcn, grad_ext_fcn)
...
...     m.obj = pyo.Objective(
...         expr=(m.z[0]-1.0)**2 + (m.z[0]-m.z[1])**2 + (m.z[2]-1.0)**2 \
...            + (m.x[0]-1.0)**4 + (m.x[1]-1.0)**6
...     )
...
...     m.c1 = pyo.Constraint(
...         expr=m.x[0] * m.z[0]**2 + m.ext_fcn(m.x[0], m.x[1]) == 2*pyo.sqrt(2.0)
...         )
...     m.c2 = pyo.Constraint(expr=m.z[2]**4 * m.z[1]**2 + m.z[1] == 8+pyo.sqrt(2.0))
...     return m
>>> model = create_model()

Step 3: Solve with TRF

Note

Reminder from earlier that the solve method requires the user pass the model and a list of variables which represent the degrees of freedom in the model. The user may also pass a low-fidelity/ideal model (or “basis function”) to this method to improve convergence.

>>> # === Instantiate the TRF solver object ===
>>> trf_solver = pyo.SolverFactory('trustregion')
>>> # === Solve with TRF ===
>>> result = trf_solver.solve(model, [model.z[0], model.z[1], model.z[2]])
EXIT: Optimal solution found.
...

The solve method returns a clone of the original model which has been run through TRF algorithm, thus leaving the original model intact.

Warning

TRF is still under a beta release. Please provide feedback and/or report any problems by opening an issue on the Pyomo GitHub page.