Revisiting the Beam Problem - Minimizing Stress with KS Constraints and BSplines

The following example shows the optimization of a simple beam with rectangular cross section. The beam has been subdivided into elements along the length of the beam, and one end is fixed. The goal is to minimize the volume (and hence the mass of the homogeneous beam) by varying the thickness in each element without exceeding a maximum stress constraint while the beam is subject to multiple load cases, each one being a distributed force load that varies sinusoidally along the span.

Constraining the bending stress on each element leads to a more computationally expensive derivative calculation, so we will use the KSFunction to reduce the stress vector for each load case to a single scalar value. To do so, we also need to insert an ExecComp component that converts the stress into a form where a negative value means it is satisfied, and a positive value means it is violated.

The problem presented here is also an example of a multi-point implementation, where we create a separate instance of the parts of the calculation that are impacted by different load cases. This enables our model to take advantage of multiple processors when run in parallel.

If we allow the optimizer to vary the thickness of each element, then we have a design variable vector that is as wide as the number of elements in the model. This may perform poorly if we have a large number of elements in the beam. If we assume that the optimal beam thickness is going to have a smooth continuous variation over the length of the beam, then it is a good candidate for using an interpolation component like SplineComp to reduce the number of design variables we need.

For this example, we have 25 elements, but can reduce that to 5 control points for the optimizer’s design variables by including the SplineComp.

The code for the top system is this:

"""
This is a multipoint implementation of the beam optimization problem.

This version minimizes volume while satisfying a max bending stress constraint in each element
for each loadcase.
"""
import numpy as np

from openmdao.test_suite.test_examples.beam_optimization.components.local_stiffness_matrix_comp import LocalStiffnessMatrixComp
from openmdao.test_suite.test_examples.beam_optimization.components.moment_comp import MomentOfInertiaComp
from openmdao.test_suite.test_examples.beam_optimization.components.multi_states_comp import MultiStatesComp
from openmdao.test_suite.test_examples.beam_optimization.components.multi_stress_comp import MultiStressComp
from openmdao.test_suite.test_examples.beam_optimization.components.volume_comp import VolumeComp
from openmdao.utils.spline_distributions import sine_distribution


def divide_cases(ncases, nprocs):
    """
    Divide up load cases among available procs.

    Parameters
    ----------
    ncases : int
        Number of load cases.
    nprocs : int
        Number of processors.

    Returns
    -------
    list of list of int
        Integer case numbers for each proc.
    """
    data = []
    for j in range(nprocs):
        data.append([])

    wrap = 0
    for j in range(ncases):
        idx = j - wrap
        if idx >= nprocs:
            idx = 0
            wrap = j

        data[idx].append(j)

    return data


class MultipointBeamGroup(om.Group):
    """
    System setup for minimization of volume (i.e., mass) subject to KS aggregated bending stress constraints.
    """

    def initialize(self):
        self.options.declare('E')
        self.options.declare('L')
        self.options.declare('b')
        self.options.declare('volume')
        self.options.declare('max_bending')
        self.options.declare('num_elements', 5)
        self.options.declare('num_cp', 50)
        self.options.declare('num_load_cases', 1)
        self.options.declare('parallel_derivs', False, types=bool, allow_none=True)
        self.options.declare('ks_add_constraint', default=False, types=bool)

    def setup(self):
        E = self.options['E']
        L = self.options['L']
        b = self.options['b']
        volume = self.options['volume']
        max_bending = self.options['max_bending']
        num_elements = self.options['num_elements']
        num_nodes = num_elements + 1
        num_cp = self.options['num_cp']
        num_load_cases = self.options['num_load_cases']
        parallel_derivs = self.options['parallel_derivs']

        x_interp = sine_distribution(num_elements)
        comp = om.SplineComp(method='bsplines', num_cp=num_cp, x_interp_val=x_interp)
        comp.add_spline(y_cp_name='h_cp', y_interp_name='h')
        self.add_subsystem('interp', comp)

        I_comp = MomentOfInertiaComp(num_elements=num_elements, b=b)
        self.add_subsystem('I_comp', I_comp)

        comp = LocalStiffnessMatrixComp(num_elements=num_elements, E=E, L=L)
        self.add_subsystem('local_stiffness_matrix_comp', comp)

        # Parallel Subsystem for load cases.
        par = self.add_subsystem('parallel', om.ParallelGroup())

        # Determine how to split cases up over the available procs.
        nprocs = self.comm.size
        divide = divide_cases(num_load_cases, nprocs)

        for j, this_proc in enumerate(divide):
            num_rhs = len(this_proc)

            name = 'sub_%d' % j
            sub = par.add_subsystem(name, om.Group())

            # Load is a sinusoidal distributed force of varying spatial frequency.
            force_vector = np.zeros((2 * num_nodes, num_rhs))
            for i, k in enumerate(this_proc):

                end = 1.5 * np.pi
                if num_load_cases > 1:
                    end += k * 0.5 * np.pi / (num_load_cases - 1)

                x = np.linspace(0, end, num_nodes)
                f = - np.sin(x)
                force_vector[0:-1:2, i] = f

            comp = MultiStatesComp(num_elements=num_elements, force_vector=force_vector,
                                   num_rhs=num_rhs)
            sub.add_subsystem('states_comp', comp)

            comp = MultiStressComp(num_elements=num_elements, E=E, num_rhs=num_rhs)
            sub.add_subsystem('stress_comp', comp)

            self.connect('local_stiffness_matrix_comp.K_local',
                         'parallel.%s.states_comp.K_local' % name)

            for k in range(num_rhs):
                sub.connect('states_comp.d_%d' % k,
                            'stress_comp.displacements_%d' % k,
                            src_indices=np.arange(2 *num_nodes))

                if parallel_derivs:
                    color = 'red_%d' % k
                else:
                    color = None

                comp = om.KSComp(width=num_elements, upper=max_bending,
                                 add_constraint=self.options['ks_add_constraint'],
                                 parallel_deriv_color=color)

                sub.add_subsystem('KS_%d' % k, comp)

                sub.connect('stress_comp.stress_%d' % k,
                            'KS_%d.g' % k)

                if not self.options['ks_add_constraint']:
                    sub.add_constraint('KS_%d.KS' % k, upper=0.0,
                                       parallel_deriv_color=color)

        comp = VolumeComp(num_elements=num_elements, b=b, L=L)
        self.add_subsystem('volume_comp', comp)

        self.connect('interp.h', 'I_comp.h')
        self.connect('interp.h', 'volume_comp.h')
        self.connect('I_comp.I', 'local_stiffness_matrix_comp.I')

        self.add_design_var('interp.h_cp', lower=1e-2, upper=10.)
        self.add_objective('volume_comp.volume')

Next we run the model, and choose ScipyOptimizeDriver SLSQP to be our optimizer. At the conclusion of optimization, we print out the design variable, which is the thickness for each element.

import numpy as np

from openmdao.test_suite.test_examples.beam_optimization.multipoint_beam_stress import MultipointBeamGroup

E = 1.
L = 1.
b = 0.1
volume = 0.01
max_bending = 100.0

num_cp = 5
num_elements = 25
num_load_cases = 2

model = MultipointBeamGroup(E=E, L=L, b=b, volume=volume, max_bending = max_bending,
                            num_elements=num_elements, num_cp=num_cp,
                            num_load_cases=num_load_cases)

prob = om.Problem(model=model)

prob.driver = om.ScipyOptimizeDriver()
prob.driver.options['optimizer'] = 'SLSQP'
prob.driver.options['tol'] = 1e-9
prob.driver.options['disp'] = True

prob.setup(mode='rev')

prob.run_driver()

print(prob['interp.h'][0])
Optimization terminated successfully.    (Exit mode 0)
            Current function value: 0.03494890159976476
            Iterations: 31
            Function evaluations: 51
            Gradient evaluations: 31
Optimization Complete
-----------------------------------
[0.45632322 0.45612552 0.45543324 0.45397058 0.45134629 0.44714397
 0.4410258  0.43283139 0.42265378 0.41087801 0.39817311 0.38543581
 0.37369202 0.36342186 0.35289066 0.34008776 0.32362887 0.30300358
 0.27867837 0.25204063 0.22519409 0.20063906 0.18088818 0.16807857
 0.16364105]