[PYTHON] Latin super square sampling with OpenMDAO

What is Latin Super Square Sampling (LHS)?

It is one of the multivariable stratified sampling methods. When n experiments are specified, each variable is divided into n sections, and the values are randomly extracted from the sections and combined randomly. .. I'm not sure if it's a sentence, so there was a module called pyDOE, so I'll try it.

Installation

pip install pyDOE

What is the experimental pattern when extracting the center of the interval with 5 experiments and 2 variables?

>from pyDOE import lhs
>lhs(2,5,"c")
array([[ 0.3,  0.7],
       [ 0.1,  0.1],
       [ 0.5,  0.9],
       [ 0.9,  0.5],
       [ 0.7,  0.3]])

If you run the above again, the combination will be different. Furthermore, if the argument " c " is not specified, it will be randomly extracted from the interval. It looks like this.

Latin super-square sampling of paraboloids using OpenMDAO

\begin{align}
& f(x,y) = (x-3)^2 + xy + (y+4)^2 - 3 \\
    {\rm subject \: to} \: \: \:&  -50.0\leq x \leq 50.0 \\
                                &  -50.0\leq y \leq 50.0
\end{align}

Component preparation

Prepare the following paraboloid.py file

paraboloid.py

from openmdao.api import Component
class Paraboloid(Component):
    """ Evaluates the equation f(x,y) = (x-3)^2 + xy + (y+4)^2 - 3 """
    def __init__(self):
        super(Paraboloid, self).__init__()
        self.add_param('x', val=0.0)
        self.add_param('y', val=0.0)
        self.add_output('f_xy', shape=1)
    def solve_nonlinear(self, params, unknowns, resids):
        """f(x,y) = (x-3)^2 + xy + (y+4)^2 - 3 """
        x = params['x']; y = params['y']
        unknowns['f_xy'] = (x-3.0)**2 + x*y + (y+4.0)**2 - 3.0

Problem settings

Prepare doe_paraboloid.py below. The difference from Solving the paraboloid minimization problem with OpenMDAO is outlined later.

doe_paraboloid.py

#!/bin/pyhton

from openmdao.api import IndepVarComp, Group, Problem,  SqliteRecorder
from paraboloid import Paraboloid
from openmdao.drivers.latinhypercube_driver import OptimizedLatinHypercubeDriver
from openmdao.core.mpi_wrap import MPI

if MPI: # pragma: no cover
    # if you called this script with 'mpirun', then use the petsc data passing
    from openmdao.core.petsc_impl import PetscImpl as impl
else:
    # if you didn't use `mpirun`, then use the numpy data passing
    from openmdao.api import BasicImpl as impl

top = Problem(impl=impl)
root = top.root = Group()

root.add('p1', IndepVarComp('x', 50.0), promotes=['x'])
root.add('p2', IndepVarComp('y', 50.0), promotes=['y'])
root.add('comp', Paraboloid(), promotes=['x', 'y', 'f_xy'])

top.driver = OptimizedLatinHypercubeDriver(num_samples=100, seed=0, population=20, \
            generations=4, norm_method=2, num_par_doe=5)
top.driver.add_desvar('x', lower=-50.0, upper=50.0)
top.driver.add_desvar('y', lower=-50.0, upper=50.0)

doe_paraboloid.py continued

top.driver.add_objective('f_xy')

recorder = SqliteRecorder('doe_paraboloid')
recorder.options['record_params'] = True
recorder.options['record_unknowns'] = True
recorder.options['record_resids'] = False
top.driver.add_recorder(recorder)

top.setup()
top.run()

top.cleanup()

When adding Paraboloid () to root, p1.x and comp.x etc. are automatically connected by using promotes as an argument.

This time, I used the Optimized Latin Hypercube Driver, which is a driver that makes it feel good to use GA to arrange the experimental points of LHS. OptimizedLatinHypercubeDriver has 100 samples, 0 random seeds, 20 GA individuals, 4 generations, a norm that normalizes the distance between each DOE point in GA (used in np.linalg.norm), parallel with num_par_doe The number of conversions is specified.

Do the following in the terminal

mpirun -np 5 python doe_paraboloid.py

Loading results

IPython

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import sqlitedict

db =sqlitedict.SqliteDict("doe_paraboloid","iterations")
res = np.array([[0.,0.,0.]] * len(db.items()))
for i, v in enumerate(db.items()):
    res[i,0] = v[1]["Unknowns"]["x"]
    res[i,1] = v[1]["Unknowns"]["y"]
    res[i,2] = v[1]["Unknowns"]["f_xy"]

x = np.arange(-50, 50, 4)
y = np.arange(-50, 50, 4)
X, Y = np.meshgrid(x, y)
Z = (X-3.0)**2 + X*Y + (Y+4.0)**2 - 3.0

fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,Z)
ax.scatter3D(res[:,0],res[:,1],res[:,2],c="r", marker="o")
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f_xy')
plt.show()