Using SciPy and Other Third-party Python Packages

This example demonstrates using an optimization routine from scipy to solve an inverse elastic problem.

A constant force, f, is applied to the tip of a cantilever beam. From a given beam geometry and displacement we use PFC to find the modulus of the beam. The beam is modeled with parallel bonded PFC particles.

The deflection is known analytically,

$\Delta y = \frac {f x^{2}}{2EI} \left(L - \frac{x}{3}\right)$

where f is the applied force, E is the Youngs Modulus, I is the moment of inertia and x is the length along the beam and L is the total length of the beam.

This is a simplification of a PFC verification problem see Cantilever Beam for more information.

import matplotlib

matplotlib.use("QT5Agg")

matplotlib.rcParams['backend.qt5']='PySide2'

import pylab as plt

import numpy as np

from scipy.constants import inch

import itasca as it

it.command("python-reset-state false")

from itasca import ballarray as ba

r = 8.0  # beam radius [inches]

l = 200.0 # beam length [inches]

x = np.linspace(0, l, 100)

E = 1e4 # [ksi]

I = 0.25*np.pi*r**4 # moment of inertia.

f = 1.0 # [kip] force applied in the positive y direction


First we import the required libraries and setup some parameters for the beam.

def analytical_deflection(x):

"Return the analytical solution for y deflection at the x point(s)"

return f * x**2 / 2.0 / E / I * ( l - x / 3.0)

def pfc_deflection(E_trial):

"""For a given parallel bond stiffness, return a tuple containing the

ball x positions and the y displacements.

calling this function cycles the PFC engine

"""

n = 11 # number of pfc particles

pb_kn = E_trial/l*(n-1)

print ("calculating deflection with pb_kn = {} ... ".format(pb_kn),)

it.command("""

model new

model large-strain on

model random

model domain ext -50 250 -50 50

contact cmat default model linearpbond

""")

for i in range(n):

it.ball.create(10.0, (i*20.0, 0.0, 0.0))

tip_ball = it.ball.find(n)

tip_ball.set_force_app_y(f)

it.command("""

ball attribute dens=1000 damp 0.7

ball property 'kn'=0 'ks'=0

model clean

contact method bond

contact property pb_ten 1e5 pb_coh 1e5 pb_kn {kn} pb_ks 188 pb_rmul 0.8

ball fix velocity-x velocity-y velocity-z ...

spin-x spin-y spin-z range id=1

model solve ratio-maximum 1e-4

""".format(kn=pb_kn))

pfc_x = ba.pos()[:,0]

pfc_dy = ba.disp()[:,1]

print ("done")

return pfc_x, pfc_dy


Two functions are defined one to give the analytical solution and one to give the PFC solution.

We can test everything is working by calling the pfc_deflection function with the known correct modulus. The pfc and analytical solution should match,

pfc_x, pfc_dy = pfc_deflection(E)

plt.plot(x, analytical_deflection(x))

plt.plot(pfc_x, pfc_dy, "o")

plt.xlabel("Length along beam [inches]")

plt.ylabel("Beam vertical displacement [inches]")

plt.draw()

plt.savefig("p3d-opt_comp.png")

output:
calculating deflection with pb_kn = 500.0 ...
done


The PFC beam deflection is shown as points. The PFC solution matches well against the analytical solution which is shown as a solid blue line.

Next we do the actual optimization calculation. The function fminbound is imported from the scipy.optimize module. This function takes three arguments, a callable function, an upper bound and a lower bound. We define the function residual to return the sum of the square of the difference between the PFC and analytical solutions. We bound the optimization between 1e2 and 5.5e4.

from scipy.optimize import fminbound

def residual(pb_kn):

"""Return the sum of squares of the residual. This is a scalar measure

of the difference between the PFC solution and the analytical solution.

"""

pfc_x, pfc_dy = pfc_deflection(pb_kn)

a_dy = analytical_deflection(pfc_x)

plt.close('all')

plt.plot(x, analytical_deflection(x))

plt.plot(pfc_x, pfc_dy, "o")

plt.xlabel("Length along beam [inches]")

plt.ylabel("Beam vertical displacement [inches]")

plt.draw()

return ((a_dy - pfc_dy)**2).sum()

E_opt =  fminbound(residual, 1e2, 5.5e4)

print ("""

E specified: {}

E calculated: {}""".format(E, E_opt))

plt.close('all')

output:
calculating deflection with pb_kn = 1053.4967008815386 ...
done
calculating deflection with pb_kn = 1701.5032991184614 ...
done
calculating deflection with pb_kn = 653.0065982369227 ...
done
calculating deflection with pb_kn = 405.490102644616 ...
done
calculating deflection with pb_kn = 530.7543627593386 ...
done
calculating deflection with pb_kn = 529.6311596321268 ...
done
calculating deflection with pb_kn = 482.21349526223526 ...
done
calculating deflection with pb_kn = 501.8190835557665 ...
done
calculating deflection with pb_kn = 501.70523214550855 ...
done
calculating deflection with pb_kn = 499.7989842280888 ...
done
calculating deflection with pb_kn = 493.08192515191905 ...
done
calculating deflection with pb_kn = 500.68789623953967 ...
done
calculating deflection with pb_kn = 500.5925746896699 ...
done
calculating deflection with pb_kn = 500.66760871604833 ...
done
calculating deflection with pb_kn = 500.6732746302266 ...
done
calculating deflection with pb_kn = 500.6776216232916 ...
done
calculating deflection with pb_kn = 500.6759566385157 ...
done
calculating deflection with pb_kn = 500.681546177477 ...
done
calculating deflection with pb_kn = 500.67693683977063 ...
done
calculating deflection with pb_kn = 500.67762921620755 ...
done
calculating deflection with pb_kn = 500.6791253622799 ...
done
calculating deflection with pb_kn = 500.678200693155 ...
done
calculating deflection with pb_kn = 500.6778475009777 ...
done
calculating deflection with pb_kn = 500.6777125935705 ...
done
calculating deflection with pb_kn = 500.6776610635263 ...
done
calculating deflection with pb_kn = 500.6776413808009 ...
done

E specified: 10000.0
E calculated: 10013.55258432415


The calculated value for E and the known value are close.

Example File

The source code for the examples here is using_scipy.py.