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:

LocalStiffnessMatrixComp class definition

class LocalStiffnessMatrixComp(om.ExplicitComponent):

    def initialize(self):
        self.options.declare('num_elements', types=int)
        self.options.declare('E')
        self.options.declare('L')

    def setup(self):
        num_elements = self.options['num_elements']
        E = self.options['E']
        L = self.options['L']

        self.add_input('I', shape=num_elements)
        self.add_output('K_local', shape=(num_elements, 4, 4))

        L0 = L / num_elements
        coeffs = np.empty((4, 4))
        coeffs[0, :] = [12, 6 * L0, -12, 6 * L0]
        coeffs[1, :] = [6 * L0, 4 * L0 ** 2, -6 * L0, 2 * L0 ** 2]
        coeffs[2, :] = [-12, -6 * L0, 12, -6 * L0]
        coeffs[3, :] = [6 * L0, 2 * L0 ** 2, -6 * L0, 4 * L0 ** 2]
        coeffs *= E / L0 ** 3

        self.mtx = np.zeros((num_elements, 4, 4, num_elements))
        for ind in range(num_elements):
            self.mtx[ind, :, :, ind] = coeffs

        self.declare_partials('K_local', 'I',
            val=self.mtx.reshape(16 * num_elements, num_elements))

    def compute(self, inputs, outputs):
        outputs['K_local'] = 0
        for ind in range(self.options['num_elements']):
            outputs['K_local'][ind, :, :] = self.mtx[ind, :, :, ind] * inputs['I'][ind]

MultiStatesComp class definition

class MultiStatesComp(om.ImplicitComponent):

    def initialize(self):
        self.options.declare('num_elements', types=int)
        self.options.declare('force_vector', types=np.ndarray)
        self.options.declare('num_rhs', types=int)

    def setup(self):
        num_elements = self.options['num_elements']
        num_nodes = num_elements + 1
        size = 2 * num_nodes + 2
        num_rhs = self.options['num_rhs']

        self.add_input('K_local', shape=(num_elements, 4, 4))
        for j in range(num_rhs):
            self.add_output('d_%d' % j, shape=size)

        cols = np.arange(16*num_elements)
        rows = np.repeat(np.arange(4), 4)
        rows = np.tile(rows, num_elements) + np.repeat(np.arange(num_elements), 16) * 2

        for j in range(num_rhs):
            disp = 'd_%d' % j

            self.declare_partials(disp, 'K_local', rows=rows, cols=cols)
            self.declare_partials(disp, disp)

    def apply_nonlinear(self, inputs, outputs, residuals):
        num_rhs = self.options['num_rhs']

        self.K = self.assemble_CSC_K(inputs)
        for j in range(num_rhs):
            force_vector = np.concatenate([self.options['force_vector'][:, j], np.zeros(2)])
            residuals['d_%d' % j] = self.K.dot(outputs['d_%d' % j]) - force_vector

    def solve_nonlinear(self, inputs, outputs):
        num_rhs = self.options['num_rhs']

        self.K = self.assemble_CSC_K(inputs)
        self.lu = splu(self.K)

        for j in range(num_rhs):

            force_vector = np.concatenate([self.options['force_vector'][:, j], np.zeros(2)])
            outputs['d_%d' % j] = self.lu.solve(force_vector)

    def linearize(self, inputs, outputs, jacobian):
        num_rhs = self.options['num_rhs']
        num_elements = self.options['num_elements']

        self.K = self.assemble_CSC_K(inputs)
        self.lu = splu(self.K)

        i_elem = np.tile(np.arange(4), 4)
        i_d = np.tile(i_elem, num_elements) + np.repeat(np.arange(num_elements), 16) * 2

        K_dense = self.K.toarray()

        for j in range(num_rhs):
            disp = 'd_%d' % j
            jacobian[disp, 'K_local'] = outputs[disp][i_d]
            jacobian[disp, disp] = K_dense

    def solve_linear(self, d_outputs, d_residuals, mode):
        num_rhs = self.options['num_rhs']

        for j in range(num_rhs):
            disp = 'd_%d' % j

            if mode == 'fwd':
                d_outputs[disp] = self.lu.solve(d_residuals[disp])
            else:
                d_residuals[disp] = self.lu.solve(d_outputs[disp])

    def assemble_CSC_K(self, inputs):
        """
        Assemble the stiffness matrix in sparse CSC format.

        Returns
        -------
        ndarray
            Stiffness matrix as dense ndarray.
        """
        num_elements = self.options['num_elements']
        num_nodes = num_elements + 1
        num_entry = num_elements * 12 + 4
        ndim = num_entry + 4

        data = np.zeros((ndim, ), dtype=inputs._get_data().dtype)
        cols = np.empty((ndim, ))
        rows = np.empty((ndim, ))

        # First element.
        data[:16] = inputs['K_local'][0, :, :].flat
        cols[:16] = np.tile(np.arange(4), 4)
        rows[:16] = np.repeat(np.arange(4), 4)

        j = 16
        for ind in range(1, num_elements):
            ind1 = 2 * ind
            K = inputs['K_local'][ind, :, :]

            # NW quadrant gets summed with previous connected element.
            data[j-6:j-4] += K[0, :2]
            data[j-2:j] += K[1, :2]

            # NE quadrant
            data[j:j+4] = K[:2, 2:].flat
            rows[j:j+4] = np.array([ind1, ind1, ind1 + 1, ind1 + 1])
            cols[j:j+4] = np.array([ind1 + 2, ind1 + 3, ind1 + 2, ind1 + 3])

            # SE and SW quadrants together
            data[j+4:j+12] = K[2:, :].flat
            rows[j+4:j+12] = np.repeat(np.arange(ind1 + 2, ind1 + 4), 4)
            cols[j+4:j+12] = np.tile(np.arange(ind1, ind1 + 4), 2)

            j += 12

        data[-4:] = 1.0
        rows[-4] = 2 * num_nodes
        rows[-3] = 2 * num_nodes + 1
        rows[-2] = 0.0
        rows[-1] = 1.0
        cols[-4] = 0.0
        cols[-3] = 1.0
        cols[-2] = 2 * num_nodes
        cols[-1] = 2 * num_nodes + 1

        n_K = 2 * num_nodes + 2
        return coo_matrix((data, (rows, cols)), shape=(n_K, n_K)).tocsc()

MomentOfInertiaComp class definition

class MomentOfInertiaComp(om.ExplicitComponent):

    def initialize(self):
        self.options.declare('num_elements', types=int)
        self.options.declare('b')

    def setup(self):
        num_elements = self.options['num_elements']

        self.add_input('h', shape=num_elements)
        self.add_output('I', shape=num_elements)

    def setup_partials(self):
        rows = cols = np.arange(self.options['num_elements'])
        self.declare_partials('I', 'h', rows=rows, cols=cols)

    def compute(self, inputs, outputs):
        outputs['I'] = 1./12. * self.options['b'] * inputs['h'] ** 3

    def compute_partials(self, inputs, partials):
        partials['I', 'h'] = 1./4. * self.options['b'] * inputs['h'] ** 2

MultiStressComp class definition

class MultiStressComp(om.ExplicitComponent):

    def initialize(self):
        self.options.declare('num_elements', types=int)
        self.options.declare('num_rhs', types=int)
        self.options.declare('E')

    def setup(self):
        num_elements = self.options['num_elements']
        num_nodes = num_elements + 1
        num_rhs = self.options['num_rhs']

        self.add_input('h', shape=num_elements)

        for j in range(num_rhs):
            self.add_input('displacements_%d' % j, shape=2 * num_nodes)
            self.add_output('stress_%d' % j, shape=num_elements)

            self.declare_partials(of='stress_%d' % j, wrt='displacements_%d' % j)
            self.declare_partials(of='stress_%d' % j, wrt='h')

    def compute(self, inputs, outputs):
        num_rhs = self.options['num_rhs']
        tk = inputs['h'] * 0.5
        E = self.options['E']

        for j in range(num_rhs):
            ang = inputs['displacements_%d' % j][1::2]
            d_ang = ang[1:] - ang[0:-1]
            outputs['stress_%d' % j] = tk * E * d_ang

    def compute_partials(self, inputs, partials):
        num_elements = self.options['num_elements']
        num_nodes = num_elements + 1
        num_states = num_nodes * 2
        num_rhs = self.options['num_rhs']
        tk = inputs['h'] * 0.5
        E = self.options['E']

        for j in range(num_rhs):
            ang = inputs['displacements_%d' % j][1::2]
            d_ang = ang[1:] - ang[0:-1]
            partials['stress_%d' % j, 'h'] = np.diag(0.5 * E * d_ang)

            J = np.zeros((num_elements, num_states))
            J[:, 1:-1:2] = -np.diag(tk * E)
            J[:, 3::2] += np.diag(tk * E)
            partials['stress_%d' % j, 'displacements_%d' % j] = J

VolumeComp class definition

class VolumeComp(om.ExplicitComponent):

    def initialize(self):
        self.options.declare('num_elements', types=int)
        self.options.declare('b', default=1.)
        self.options.declare('L')

    def setup(self):
        num_elements = self.options['num_elements']
        b = self.options['b']
        L = self.options['L']
        L0 = L / num_elements

        self.add_input('h', shape=num_elements)
        self.add_output('volume')

        self.declare_partials('volume', 'h', val=b * L0)

    def compute(self, inputs, outputs):
        L0 = self.options['L'] / self.options['num_elements']

        outputs['volume'] = np.sum(inputs['h'] * self.options['b'] * L0)

"""
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
import openmdao.api as om

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_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.

from openmdao.test_suite.test_examples.beam_optimization.components.multi_states_comp import MultiStatesComp

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.03494890159977874
            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]