Optimizing the Thickness Distribution of a Cantilever Beam Using the Adjoint Method

In this example, we optimize the thickness (height) distribution of a cantilever beam using the adjoint method to compute the gradient. We use Euler–Bernoulli beam theory and assume a rectangular section.

Background

The optimization problem is:

\[\begin{split} \begin{array}{r c l} \text{minimize} & & f^T d \\ \text{with respect to} & & h \\ \text{subject to} & & \text{sum}(h) b L_0 = \text{volume} \\ \end{array} \end{split}\]

where \(f\) is the vector of forces, \(h\) is the vector of beam heights, and \(L_0\) is the length of a single beam element.

The displacements vector \(d\) is given by

\[ K d = f \]

where \(K\) is the stiffness matrix. However, in practice, we augment the linear system with Lagrange multipliers to apply the boundary constraints at the first node.

Since our model contains a system of equations, we use the adjoint method to compute the gradient of the objective with respect to the beam height vector. The model is shown below.

Optimization terminated successfully.    (Exit mode 0)
            Current function value: 23762.15367744602
            Iterations: 146
            Function evaluations: 401
            Gradient evaluations: 146
Optimization Complete
-----------------------------------

Implementation: list of components

There are 5 components that compute:

  1. Moment of inertia for each element

  2. Local stiffness matrix for each element

  3. Displacements from solution of the \(Kd=f\) linear system augmented with the Lagrange multipliers

  4. Compliance

  5. Volume

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
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 = 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]
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import splu

class StatesComp(om.ImplicitComponent):

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

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

        self.add_input('K_local', shape=(num_elements, 4, 4))
        self.add_output('d', 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

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

    def apply_nonlinear(self, inputs, outputs, residuals):
        force_vector = np.concatenate([self.options['force_vector'], np.zeros(2)])

        self.K = self.assemble_CSC_K(inputs)
        residuals['d'] = self.K.dot(outputs['d'])  - force_vector

    def solve_nonlinear(self, inputs, outputs):
        force_vector = np.concatenate([self.options['force_vector'], np.zeros(2)])

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

        outputs['d'] = self.lu.solve(force_vector)

    def linearize(self, inputs, outputs, jacobian):
        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

        jacobian['d', 'K_local'] = outputs['d'][i_d]

        jacobian['d', 'd'] = self.K.toarray()

    def solve_linear(self, d_outputs, d_residuals, mode):
        if mode == 'fwd':
            d_outputs['d'] = self.lu.solve(d_residuals['d'])
        else:
            d_residuals['d'] = self.lu.solve(d_outputs['d'])

    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()
class ComplianceComp(om.ExplicitComponent):

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

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

        self.add_input('displacements', shape=2 * num_nodes)
        self.add_output('compliance')

    def setup_partials(self):
        num_nodes = self.options['num_elements'] + 1
        force_vector = self.options['force_vector']
        self.declare_partials('compliance', 'displacements',
                              val=force_vector.reshape((1, 2 * num_nodes)))

    def compute(self, inputs, outputs):
        outputs['compliance'] = np.dot(self.options['force_vector'], inputs['displacements'])
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)

Implementation: Optimization Script

Here is the optimization script:

from openmdao.test_suite.test_examples.beam_optimization.beam_group import BeamGroup

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

num_elements = 50

prob = om.Problem(model=BeamGroup(E=E, L=L, b=b, volume=volume, num_elements=num_elements))

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

prob.setup()

prob.run_driver()

print(prob['h'])
Optimization terminated successfully.    (Exit mode 0)
            Current function value: 23762.15367744602
            Iterations: 146
            Function evaluations: 401
            Gradient evaluations: 146
Optimization Complete
-----------------------------------
[0.14915751 0.14764323 0.14611327 0.14456711 0.14300423 0.14142409
 0.13982612 0.13820964 0.13657404 0.13491862 0.13324263 0.1315453
 0.12982578 0.12808319 0.12631654 0.12452485 0.122707   0.12086182
 0.11898803 0.11708426 0.11514903 0.11318071 0.11117757 0.10913767
 0.10705892 0.10493901 0.1027754  0.10056527 0.09830547 0.0959925
 0.09362244 0.09119083 0.0886926  0.08612199 0.0834723  0.08073574
 0.07790321 0.07496383 0.07190452 0.0687093  0.06535832 0.06182633
 0.05808048 0.05407661 0.04975295 0.0450185  0.03972912 0.03363155
 0.02620192 0.01610863]

optimized