How to best setup OpenMDAO for a multi-fidelity optimisation problem

Question

I am setting up a multi-fidelity optimiztion routine using OpenMDAO, with all computations done by a high-fidelity solver and gradient computations done by a low-fidelity solver.

So far I have set it up with three components in a group: component 1 performs the high-fidelity calculations; component 2 performs the low-fidelity calculations at three design points - one point is used for calculation of a tunning factor between high and low fidelity, and the other two are used to compute the derivatives with a central difference (coded "by hand"). The outputs of component 1 (some absolute responses) and derivatives from component 2 are then passed to component 3, which defines design objectives and constraints. In this component, I declare the partials of objectives and constraints w.r.t. design variables as analytic, but then in the compute_partials I simply set them to the values passed from component 2.

Is there a more elegant way to handle this kind of setup in OpenMDAO?

I would like to also ensure that OpenMDAO does not compute other derivatives and only uses the ones provided "analytically". My understanding is that if I skip the setup_partials on components 1 and 2, these will not be computed in any way. Is that correct?

Happy to provide more information if needed.

I would appreciate any advice.

Answer 1

I don't think skipping the setup_partials on the two upstream components will work here. OpenMDAO's derivatives solvers, in the absence of any solver loops or cyclic data connections, can be thought of as equivalent to the chain rule. By leaving out the setup_partials of the upstream components you're essentially setting those to 0, so the chain rule will zero everything out.

OpenMDAO wasn't designed to support the approach you're using here. Its definitely possible to trick it into doing it by setting some partials in the upstream components to 1 (so the chain rule basically ignores them), but the result would be a dirty hack and I hesitate to even write up a sample. You would be able to get "correct" totals, but all your partials would look wrong and I just think its a bad idea.

Instead, I'd suggest you combine your high/low fidelity solvers into a single OM component. Implement both the high and low solves in compute and then do the necessary derivative calculations in the compute_partials (or implicit component equivalents).

If you provide a sample toy script, Id be willing to modify it to show you what I mean in more detail. But generally, "hidding" the multi-fidelity behavior all inside a single component is the better approach IMO.

Answer 2

Based on the response, I have modified the approach to reflect my understanding of how it can be improved, please see the toy problem below (exceeds the length limit for an edit).

It would be better to let OpenMDAO handle the derivatives in some way, instead of the custom calculation, but I have trouble figuring out how I could make it compute the derivatives of objective/constraints wrt design variables based on LF solution and the objective/constraints themselves based on a HF solution? I would appreciate further advice.

import os
import numpy as np
import openmdao.api as om



step = 0.01

def finite_diff(v, response, step):
    J = np.zeros(len(v))
    for i in range(len(v)):
        partial = (response['plus_delta'][i] - response['minus_delta'][i]) / 2*step
        J[i] = partial
        
    return J

class Main(om.ExplicitComponent):
    """All in one component"""
    
    def setup(self):
        # inputs
        self.add_input('v', shape=(3), val=np.array([100,10,8]), units="m", desc="Design variables")
        
        # Outputs
        self.add_output('R1_HF', 0.0, units="m", desc="R1 obtained by HF solve")
        self.add_output('R2_HF', 0.0, units="m", desc="R2 obtained by HF solve")    
        
    def setup_partials(self):
        self.declare_partials('R1_HF', 'v', method='exact')
        self.declare_partials('R2_HF', 'v', method='exact')
        
    def compute(self, inputs, outputs):
         
        # define inputs
        v = inputs['v']
        
        #  HF computations
        global R1_HF
        global R2_HF
        
        R1_HF = .9*(v[0]+v[2])
        R2_HF = .9*(v[1]-v[2])
        
        # define outputs
        outputs['R1_HF']   = R1_HF
        outputs['R2_HF']   = R2_HF 
        
    def compute_partials(self, inputs, J):
        """Jacobian of partial derivatives"""
        
        # define inputs
        v = inputs['v']
        
        #  LF computations
        R1_LF = {}
        R2_LF = {}
        
        # center point
        R1_LF['middle'] = v[0]+v[2]
        R2_LF['middle'] = v[1]-v[2]
                
        # tunning factor(s)
        f1 = R1_HF/R1_LF['middle']
        f2 = R2_HF/R2_LF['middle']
        
        # perturbed design variables (2N cases for N vars)        
        R1_LF['minus_delta'] = []
        R1_LF['plus_delta'] = []
        R2_LF['minus_delta'] = []
        R2_LF['plus_delta'] = []
        for i in range(len(v)):
            
            # zero-out delta vector for next point
            delta = [0,0,0]   
            
            # positive perturbation
            delta[i] = step
            dv1 = v + delta
            
            # plus delta point
            R1_LF['plus_delta'].append(v[0]+v[2]-.1)
            R2_LF['plus_delta'].append(v[0]-v[2]-.11)
        
            # negative perturbation
            dv2 = v - delta
            
            # plus delta point
            R1_LF['minus_delta'].append(v[0]+v[2]-.15)       
            R2_LF['minus_delta'].append(v[0]-v[2]-.16)       

        # tune the response(s) based on HF results
        for key in R1_LF.keys():
            if not key == 'middle':
                R1_LF[key] = list(np.array(R1_LF[key]) * f1)
                R2_LF[key] = list(np.array(R2_LF[key]) * f2)
            else:
                R1_LF[key] = np.array(R1_LF[key]) * f1
                R2_LF[key] = np.array(R2_LF[key]) * f2
        
        # define outputs           
        J['R1_HF', 'v']  = finite_diff(v, R1_LF, step)
        J['R2_HF', 'v'] = finite_diff(v, R2_LF, step)

class Group(om.Group):
    def setup(self):
        self.add_subsystem('Main_system', Main(), promotes_inputs=['v'],
                    promotes_outputs=['R1_HF','R2_HF']) 


if __name__ == "__main__":

    # initialize
    prob = om.Problem(model=Group())
    model = prob.model

    # setup the model
    model.add_design_var('v', lower=[80, 4, 4], upper=[150, 16, 16])
    model.add_objective('R1_HF')
    model.add_constraint('R2_HF', lower=-1e9, upper=1e9)

    # >>> check partials, for debugging
    prob.setup(mode='rev')
    prob.run_model()
    data = prob.check_partials(compact_print=True, includes='Main_system')