KSComp with variable-size constraint vector input

Question

Summary

I'm trying to perform structural optimization coupling Nastran SOL 106 and OpenMDAO. I want to minimize mass subject to nonlinear stress and stability constraints. My structure is a box beam reinforced with ribs and stiffeners clamped at the root and loaded with a concetrated force at the tip. At the moment I am only using a single design variable, that is the wall thickness (same for all elements of the structure).

For both the stress and the stability constraint, I would like to use OpenMDAO's KSComp component. However, I have a problem with the fixed width of the constraint vector input. While for the stress aggregation function this is not a problem, as the number of elements is fixed and I can easily calculate it, the vector of the stability constraint may change size at every iteration. This is because the constraint is imposed on the N lowest eigenvalues of the tangent stiffness matrix for each converged iteration of the nonlinear analysis (they have to be larger than 0), where the nonlinear solver is an adaptive arc-length method. This means that the number of iterations is not known a priori and that it usually changes for every analysis. As a consequence, I have a N x no_iterations array that I want to flatten and feed to the KSComp component, but I can't do that because the number of iterations is variable and cannot be fixed.

What I've tried

I have created an explicit component and a group as you see below. Inside my explicit component Sol106Comp I set the value of thickness, run the Nastran analysis, calculate the functions of interest and then call the method compute_ks_function that I have defined taking as a reference the source code of KSComp. This works, but I would like to use KSComp to benefit of its membership to the OpenMDAO world.

From this answer, it looks like KSComp is not meant at all to have a dynamic width of the contraint vector input. Is this still the case?

I've read that a solution may be over-allocation, and it may be even natural in my case, as I know the maximum number of iterations in the nonlinear analysis, so I know the maximum possible size of the constraint vector input. However what value should I give to the inactive entries of the constraint vector? In my case the constraint is satisfied when the elements of the constraint vectors are larger then 0, so I cannot set the inactive entries to 0. I could use a very large value, but wouldn't this affect my KS function in some way?

Code

import openmdao.api as om
import numpy as np
from resources import box_beam_utils, pynastran_utils
import os
from pyNastran.bdf.mesh_utils.mass_properties import mass_properties
from pyNastran.bdf.bdf import BDF


class Sol106Comp(om.ExplicitComponent):
    """
    Evaluates the mass of the model, the von mises stresses and the lowest eigenvalues of the tangent stiffness matrix.
    """

    def initialize(self):
        self.options.declare('box_beam_bdf', types=BDF)
        self.options.declare('sigma_y', types=float)

    def setup(self):
        self.add_input('t')
        self.add_output('mass')
        self.add_output('ks_stress')
        self.add_output('ks_stability')
        self.add_discrete_output('sol_106_op2', None)
        self.add_discrete_output('eigenvalues', None)

    def setup_partials(self):
        # Finite difference all partials
        self.declare_partials('*', '*', method='fd')

    def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
        """
        Run SOL 106 and calculate output functions
        """
        # Assign variables
        box_beam_bdf = self.options['box_beam_bdf']
        yield_strength = self.options['sigma_y']
        # Assign thickness to PSHELL card
        box_beam_bdf.properties[1].t = inputs['t'][0]
        # Set up directory and input name
        analysis_directory_name = 'Optimization'
        analysis_directory_path = os.path.join(os.getcwd(), 'analyses', analysis_directory_name)
        input_name = 'box_beam_sol_106'
        # Run Nastran and return OP2 file
        sol_106_op2 = pynastran_utils.run_tangent_stiffness_matrix_eigenvalue_calculation(bdf_object=box_beam_bdf.__deepcopy__({}), method_set_id=21, no_eigenvalues=10,
                                                                                         analysis_directory_path=analysis_directory_path, input_name=input_name, run_flag=True)
        # Calculate mass
        outputs['mass'] = mass_properties(box_beam_bdf)[0]
        # Find von mises stresses and aggregate with KS function
        stresses = sol_106_op2.nonlinear_cquad4_stress[1].data[-1, :, 5]
        outputs['ks_stress'] = self.compute_ks_function(stresses, upper=yield_strength)
        # Read eigenvalues of tangent stiffness matrix and aggregate with KS function
        f06_filepath = os.path.join(analysis_directory_path, input_name + '.f06')  # path to .f06 file
        eigenvalues = pynastran_utils.read_kllrh_lowest_eigenvalues_from_f06(f06_filepath)  # read eigenvalues of kllrh matrix from f06 files
        outputs['ks_stability'] = self.compute_ks_function(eigenvalues[~np.isnan(eigenvalues)].flatten(), lower_flag=True)
        # Save OP2 object
        discrete_outputs['sol_106_op2'] = sol_106_op2
        discrete_outputs['eigenvalues'] = eigenvalues
    
    @staticmethod
    def compute_ks_function(g, rho=50, upper=0, lower_flag=False):
        """
        Compute the value of the KS function for the given array of constraints.
        Reference: https://openmdao.org/newdocs/versions/latest/_modules/openmdao/components/ks_comp.html
        """
        con_val = g - upper
        if lower_flag:
            con_val = -con_val
        g_max = np.max(np.atleast_2d(con_val), axis=-1)[:, np.newaxis]
        g_diff = con_val - g_max
        exponents = np.exp(rho * g_diff)
        summation = np.sum(exponents, axis=-1)[:, np.newaxis]
        KS = g_max + 1.0 / rho * np.log(summation)
        return KS


class BoxBeamGroup(om.Group):
    """
    System setup for minimization of mass subject to KS aggregated von mises stress and structural stability constraints.
    """

    def initialize(self):
        # Define object input variable
        ...

    def setup(self):
        # Assign variables
        ...
        # Genearte mesh
        ...
        # Create BDF object
        ...
        # Apply load
        ...
        # Create subcase
        ...
        # Set up SOL 106 analysis
        ...
        # Define SOL 106 component
        comp = Sol106Comp(box_beam_bdf=box_beam_bdf, sigma_y=yield_strength)
        self.add_subsystem('sol_106', comp)

Answer 1

You are right that overallocation is the approach you would need to take. In general, I don't like formulations where the number of constraints is changing. When that happens it means something non-differentiable has happened in your code and generally it will cause noise.

However, I understand that NASTRAN has some limitations you have to work around so either you add the KS function inside your component where you can easily handle the variable sized vector or you use overallocation. I don't think its necessary to use the OpenMDAO component if you don't have to, and in fact it may even be computationally better to skip it. You're giving OpenMDAO less intermediate variables to deal with, which generally would improve performance.

You can, as you suggest, use a very large value for the constraint. As long as its nowhere here the transition level it shouldn't cause massive numerical difficulties.

To use overallocation, you'd really have to hack the OpenMDAO KS component anyway. You'd want to put the length (i.e. number of entries to consider) into the first element of your array. Then hack the KS component to read that value (and discard it from the KS aggregation). Then use that information to pull the relevant values to aggregate and handle the derivatives appropriately. This is essentially the same thing you'd do if you just put the KS inside your component anyway.