How to optimize with two implicit, and possibly coupled, systems in OpenMDAO?

Question

I've been referring to the Sellar MDA and BalanceComp tutorials, but I'm unable to formulate an OpenMDAO architecture to optimize a coupled system of implicit equations of the type:

Minimize the function F with design variables x, y such that the following nonlinear equations are simultaneously satisfied as constraints:

R(x, y, z) = f(x, y, c) - g(x, y, d) = 0

S(x, y, z) = f(x, y, c) - h(x, y, d) = e

with constants c, d, e.

So my idea is to let F, f, g, h be ExplicitComponent classes, and let R, S be Group classes that use BalanceComp (one inside each!) to determine their respective solutions, because the systems may be used independently in other applications. The constants c, d, e can be fed via an IndepVarComp, I believe.

Solving the system together involves creating a Group component (call it, say, ImpMDA) with R, S inside it using a Group like cycle from the Sellar tutorial.

My main issue is how to "synchronize" x, y globally across the two systems in the cycle. Should R's results be fed to S's and vice versa? As they're being solved simultaneously, I would assume they're coupled.

Please note that x, y may come from some other, separate components, and hence feeding them as a vector is difficult unless some ExplicitComponent like Lister is created, which takes inputs and generates a vector to feed to ImpMDA.

While a specific implementation on a simpler toy case (if possible) would be very helpful, I would really appreciate it if I receive some general insight into modeling this problem in OpenMDAO anyway.

For the specific use-case, consider the minimization of aircraft fuel burn (objective function) with respect to variation of wing area (design variable) while trying to maintain pitch trim (implicit system 1) and pitch stability (implicit system 2) at all operating points in flight, hence determining the horizontal tail area and wing location. The center of pressure (dependent on wing area) matching the center of gravity is the requirement for pitch trim, and similarly for the neutral point for pitch stability, hence the design variable plays a role in the systems.

EDIT:

Here is a Python example of the modular approach I'm looking for:

import numpy as np
from scipy.optimize import fsolve, minimize

def f(xs, c): return -xs[0]**2 + xs[1] + c          # ExplicitComponent 1
def g(xs, d): return -xs[0]**3 + xs[1]**2 - d       # ExplicitComponent 2
def h(xs, d): return xs[0] + xs[1]**2 + d           # ExplicitComponent 3

# This function is to quickly generate residual functions shown in the question.
# The OpenMDAO implementations for different residual functions will be done either 
# via BalanceComp or ImplicitComponent, and may not share the same equation "setup" as here.
def residualer(f1, f2, a=6, b=4, const=0):           
    return lambda xs: f1(xs, a) - f2(xs, b) + const 

def F(xs):                                          # Objective function
    return xs[0]**2 - xs[1]**2 - 5

x0 = [-2, 2]                                        # Initial guess

# Problem 1
x1_opt = minimize(F, x0, 
                 method='trust-constr',
                 constraints={ 'type': 'eq', 'fun': residualer(f, g, a=1, b=2) }, 
                 options = { 'disp': False } )

print(x1_opt, '\n\n')


# Problem 2
m, n = 5, 6
def coupled(xs):                                    # Coupling the equations 
    return [residualer(f, g)(xs), residualer(f, h, a=4, b=3, const=m*n)(xs)] 

x2_opt = minimize(F, x0, 
                 method='trust-constr',
                 constraints={ 'type': 'eq', 'fun': coupled }, 
                 options = { 'disp': False } )

print(x2_opt)

which returns:

         cg_niter: 22
     cg_stop_cond: 2
           constr: [array([-8.8817842e-16])]
      constr_nfev: [102]
      constr_nhev: [0]
      constr_njev: [0]
   constr_penalty: 3.2410415627552283
 constr_violation: 8.881784197001252e-16
   execution_time: 0.028232574462890625
              fun: -10.302775637731996
             grad: array([ 1.19209290e-07, -4.60555129e+00])
              jac: [array([[-5.96046448e-08, -3.60555129e+00]])]
  lagrangian_grad: array([ 1.95345288e-07, -8.88178420e-16])
          message: '`xtol` termination condition is satisfied.'
           method: 'equality_constrained_sqp'
             nfev: 102
             nhev: 0
              nit: 23
            niter: 23
             njev: 0
       optimality: 1.9534528835387254e-07
           status: 2
          success: True
        tr_radius: 5.238777835236851e-09
                v: [array([-1.2773501])]
                x: array([1.17707733e-08, 2.30277564e+00]) 


         cg_niter: 0
     cg_stop_cond: 1
           constr: [array([4.4408921e-14, 0.0000000e+00])]
      constr_nfev: [36]
      constr_nhev: [0]
      constr_njev: [0]
   constr_penalty: 2.1867764273118655
 constr_violation: 4.440892098500626e-14
   execution_time: 0.012656688690185547
              fun: -24.59492687150548
             grad: array([  5.27636947, -10.30629829])
              jac: [array([[15.60368701, -9.30629829],
       [-6.27636948, -9.30629823]])]
  lagrangian_grad: array([1.77635684e-15, 0.00000000e+00])
          message: '`gtol` termination condition is satisfied.'
           method: 'equality_constrained_sqp'
             nfev: 36
             nhev: 0
              nit: 12
            niter: 12
             njev: 0
       optimality: 1.7763568394002505e-15
           status: 1
          success: True
        tr_radius: 4.9172480000000025
                v: [array([-0.55882674, -0.54862737])]
                x: array([2.63818474, 5.1531491 ])

Answer 1

I took a stab at doing this in a modular way. Since the script is fairly large, its on a github gist here.

Here Residualer is a group that serves the purpose of your residualer function. In its initialize method, I defined options that are specific for each instantiation.

Each instantiation creates its own subsystems for f, g, and h. In practice, I'd probably make a single component that outputs f, g, and h since they're simple equations and independent of one another.

problem1_optimizer() and problem2_optimizer() use scipy optimizer with equality constraints on the residuals, as your examples do.

The N2 for the problem2_optimizer is:

problem2_solver() uses the use_solver option on residualer to include a balance component and a NewtonSolver. This balance component adds an implicit variable for each element in x, which has a corresponding residual (R or S). These two scalar values (x1 and x2) are muxed into the x variable and then passed to the appropriate components by promotion. As coded, this requires that we have a single residual for each scalar implicit variable, and so it only works on your second example. Note the unconnected inputs in the following are due to the fact that our balance comp is only taking in one variable and using an assumed value of 0.0 for the other side of the equation. Theres no optimization in this second case since there are only two variables and two unknowns, and therefore no more degrees of freedom for optimization.