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

# OpenMDAO Partition Tree and N2 diagram.

## Implementation: group¶

We first show the Group that contains all the Component instances for the model.

from __future__ import division
import numpy as np

from openmdao.api import Group, IndepVarComp

from openmdao.test_suite.test_examples.beam_optimization.components.moment_comp import MomentOfInertiaComp
from openmdao.test_suite.test_examples.beam_optimization.components.local_stiffness_matrix_comp import LocalStiffnessMatrixComp
from openmdao.test_suite.test_examples.beam_optimization.components.global_stiffness_matrix_comp import GlobalStiffnessMatrixComp
from openmdao.test_suite.test_examples.beam_optimization.components.states_comp import StatesComp
from openmdao.test_suite.test_examples.beam_optimization.components.displacements_comp import DisplacementsComp
from openmdao.test_suite.test_examples.beam_optimization.components.compliance_comp import ComplianceComp
from openmdao.test_suite.test_examples.beam_optimization.components.volume_comp import VolumeComp

class BeamGroup(Group):

def initialize(self):

def setup(self):
num_nodes = num_elements + 1

force_vector = np.zeros(2 * num_nodes)
force_vector[-2] = -1.

inputs_comp = IndepVarComp()

I_comp = MomentOfInertiaComp(num_elements=num_elements, b=b)

comp = LocalStiffnessMatrixComp(num_elements=num_elements, E=E, L=L)

comp = GlobalStiffnessMatrixComp(num_elements=num_elements)

comp = StatesComp(num_elements=num_elements, force_vector=force_vector)

comp = DisplacementsComp(num_elements=num_elements)

comp = ComplianceComp(num_elements=num_elements, force_vector=force_vector)

comp = VolumeComp(num_elements=num_elements, b=b, L=L)

self.connect('inputs_comp.h', 'I_comp.h')
self.connect('I_comp.I', 'local_stiffness_matrix_comp.I')
self.connect(
'local_stiffness_matrix_comp.K_local',
'global_stiffness_matrix_comp.K_local')
self.connect(
'global_stiffness_matrix_comp.K',
'states_comp.K')
self.connect(
'states_comp.d',
'displacements_comp.d')
self.connect(
'displacements_comp.displacements',
'compliance_comp.displacements')
self.connect(
'inputs_comp.h',
'volume_comp.h')



## Implementation: list of components¶

There are 7 components that compute:

1. moment of inertia for each element
2. local stiffness matrix for each element
3. global stiffness matrix
4. solution of the $$Kd=f$$ linear system augmented with the Lagrange multipliers
5. extraction of just the displacements in the $$d$$ vector
6. compliance
7. volume
from __future__ import division
import numpy as np
from openmdao.api import ExplicitComponent

class MomentOfInertiaComp(ExplicitComponent):

def initialize(self):

def setup(self):

rows = np.arange(num_elements)
cols = np.arange(num_elements)
self.declare_partials('I', 'h', rows=rows, cols=cols)

def compute(self, inputs, outputs):

outputs['I'] = 1./12. * b * inputs['h'] ** 3

def compute_partials(self, inputs, partials):

partials['I', 'h'] = 1./4. * b * inputs['h'] ** 2

from __future__ import division
import numpy as np
from openmdao.api import ExplicitComponent

class LocalStiffnessMatrixComp(ExplicitComponent):

def initialize(self):

def setup(self):

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(num_elements):
outputs['K_local'][ind, :, :] = self.mtx[ind, :, :, ind] * inputs['I'][ind]

from __future__ import division
import numpy as np
from openmdao.api import ExplicitComponent

class GlobalStiffnessMatrixComp(ExplicitComponent):

def initialize(self):

def setup(self):
num_nodes = num_elements + 1

self.add_output('K', shape=(2 * num_nodes + 2, 2 * num_nodes + 2))

rows = np.zeros(16 * num_elements, int)
indices = np.arange(
((2 * num_nodes + 2) * (2 * num_nodes + 2))
).reshape((2 * num_nodes + 2, 2 * num_nodes + 2))
ind1, ind2 = 0, 0
for ind in range(num_elements):
ind2 += 16
ind1_ = 2 * ind
ind2_ = 2 * ind + 4
rows[ind1:ind2] = indices[ind1_:ind2_, ind1_:ind2_].flatten()
ind1 += 16
cols = np.arange(16 * num_elements)
self.declare_partials('K', 'K_local', val=1., rows=rows, cols=cols)

self.set_check_partial_options('K_local', step=1e0)

def compute(self, inputs, outputs):
num_nodes = num_elements + 1

outputs['K'][:, :] = 0.
for ind in range(num_elements):
ind1_ = 2 * ind
ind2_ = 2 * ind + 4

outputs['K'][ind1_:ind2_, ind1_:ind2_] += inputs['K_local'][ind, :, :]

outputs['K'][2 * num_nodes + 0, 0] = 1.0
outputs['K'][2 * num_nodes + 1, 1] = 1.0
outputs['K'][0, 2 * num_nodes + 0] = 1.0
outputs['K'][1, 2 * num_nodes + 1] = 1.0

from __future__ import division
import numpy as np
from scipy.linalg import lu_factor, lu_solve
from openmdao.api import ImplicitComponent

class StatesComp(ImplicitComponent):

def initialize(self):

def setup(self):
num_nodes = num_elements + 1
size = 2 * num_nodes + 2

rows = np.outer(np.arange(size), np.ones(size, int)).flatten()
cols = np.arange(size ** 2)
self.declare_partials('d', 'K', rows=rows, cols=cols)

self.declare_partials('d', 'd')

def apply_nonlinear(self, inputs, outputs, residuals):

residuals['d'] = np.dot(inputs['K'], outputs['d']) - force_vector

def solve_nonlinear(self, inputs, outputs):

self.lu = lu_factor(inputs['K'])

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

def linearize(self, inputs, outputs, partials):
num_nodes = num_elements + 1
size = 2 * num_nodes + 2

self.lu = lu_factor(inputs['K'])

partials['d', 'K'] = np.outer(np.ones(size), outputs['d']).flatten()
partials['d', 'd'] = inputs['K']

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

from __future__ import division
import numpy as np
from openmdao.api import ExplicitComponent

class DisplacementsComp(ExplicitComponent):

def initialize(self):

def setup(self):
num_nodes = num_elements + 1
size = 2 * num_nodes + 2

arange = np.arange(2 * num_nodes)
self.declare_partials('displacements', 'd', val=1., rows=arange, cols=arange)

def compute(self, inputs, outputs):
num_nodes = num_elements + 1

outputs['displacements'] = inputs['d'][:2 * num_nodes]

from __future__ import division
import numpy as np
from openmdao.api import ExplicitComponent

class ComplianceComp(ExplicitComponent):

def initialize(self):

def setup(self):
num_nodes = num_elements + 1

self.declare_partials('compliance', 'displacements',
val=force_vector.reshape((1, 2 * num_nodes)))

def compute(self, inputs, outputs):

outputs['compliance'] = np.dot(force_vector, inputs['displacements'])

from __future__ import division
import numpy as np
from openmdao.api import ExplicitComponent

class VolumeComp(ExplicitComponent):

def initialize(self):

def setup(self):
L0 = L / num_elements

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

def compute(self, inputs, outputs):
L0 = L / num_elements

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


## Implementation: optimization script¶

Here is the optimization script:

import numpy as np

from openmdao.api import Problem, ScipyOptimizeDriver

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 = Problem(model=BeamGroup(E=E, L=L, b=b, volume=volume, num_elements=num_elements))

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

prob.setup()

prob.run_driver()
print(prob['inputs_comp.h'])

Optimization terminated successfully.    (Exit mode 0)
Current function value: 23762.1536775
Iterations: 147
Function evaluations: 440
0.02620193  0.01610863]