Basic Tutorials
This is a collection of all tutorials demonstrating basic ParMOO functionality (collected from throughout the ParMOO User Guide).
Quickstart Demo
This is a basic example (see quickstart.py) of how to build and solve a MOOP with ParMOO, taken from the Quickstart guide.
import numpy as np
import pandas as pd
from parmoo import MOOP
from parmoo.searches import LatinHypercube
from parmoo.surrogates import GaussRBF
from parmoo.acquisitions import RandomConstraint
from parmoo.optimizers import LocalGPS
# Fix the random seed for reproducibility
np.random.seed(0)
my_moop = MOOP(LocalGPS)
my_moop.addDesign({'name': "x1",
'des_type': "continuous",
'lb': 0.0, 'ub': 1.0})
my_moop.addDesign({'name': "x2", 'des_type': "categorical",
'levels': ["good", "bad"]})
def sim_func(x):
if x["x2"] == "good":
return np.array([(x["x1"] - 0.2) ** 2, (x["x1"] - 0.8) ** 2])
else:
return np.array([99.9, 99.9])
my_moop.addSimulation({'name': "MySim",
'm': 2,
'sim_func': sim_func,
'search': LatinHypercube,
'surrogate': GaussRBF,
'hyperparams': {'search_budget': 20}})
def f1(x, s): return s["MySim"][0]
def f2(x, s): return s["MySim"][1]
my_moop.addObjective({'name': "f1", 'obj_func': f1})
my_moop.addObjective({'name': "f2", 'obj_func': f2})
def c1(x, s): return 0.1 - x["x1"]
my_moop.addConstraint({'name': "c1", 'constraint': c1})
for i in range(3):
my_moop.addAcquisition({'acquisition': RandomConstraint,
'hyperparams': {}})
my_moop.solve(5)
results = my_moop.getPF(format="pandas")
# Display solution
print(results)
# Plot results -- must have extra viz dependencies installed
from parmoo.viz import scatter
# The optional arg `output` exports directly to jpg instead of interactive mode
scatter(my_moop, output="jpeg")
The above code saves all (approximate) Pareto optimal solutions in the
results dataframe, and prints the results dataframe to the standard
output:
x1 x2 f1 f2 c1
0 0.742840 good 0.294675 0.003267 -0.642840
1 0.726092 good 0.276773 0.005462 -0.626092
2 0.605914 good 0.164766 0.037669 -0.505914
3 0.548931 good 0.121753 0.063036 -0.448931
4 0.543499 good 0.117991 0.065793 -0.443499
5 0.401011 good 0.040405 0.159192 -0.301011
6 0.353552 good 0.023578 0.199316 -0.253552
7 0.328402 good 0.016487 0.222404 -0.228402
8 0.269175 good 0.004785 0.281775 -0.169175
9 0.248183 good 0.002322 0.304502 -0.148183
And produces the following figure of the Pareto points:
Named Output Types
The following named_var_ex.py code demonstrates ParMOO’s output datatype when the
MOOP object is defined using named variables.
import numpy as np
from parmoo import MOOP
from parmoo.searches import LatinHypercube
from parmoo.surrogates import GaussRBF
from parmoo.acquisitions import UniformWeights
from parmoo.optimizers import LocalGPS
my_moop = MOOP(LocalGPS)
# Define a simulation to use below
def sim_func(x):
if x["MyCat"] == 0:
return np.array([(x["MyDes"]) ** 2, (x["MyDes"] - 1.0) ** 2])
else:
return np.array([99.9, 99.9])
# Add a design variable, simulation, objective, and constraint.
# Note the 'name' keys for each
my_moop.addDesign({'name': "MyDes",
'des_type': "continuous",
'lb': 0.0, 'ub': 1.0})
my_moop.addDesign({'name': "MyCat",
'des_type': "categorical",
'levels': 2})
my_moop.addSimulation({'name': "MySim",
'm': 2,
'sim_func': sim_func,
'search': LatinHypercube,
'surrogate': GaussRBF,
'hyperparams': {'search_budget': 20}})
my_moop.addObjective({'name': "MyObj",
'obj_func': lambda x, s: sum(s["MySim"])})
my_moop.addConstraint({'name': "MyCon",
'constraint': lambda x, s: 0.1 - x["MyDes"]})
# Extract numpy dtypes for all of this MOOP's inputs/outputs
des_dtype = my_moop.getDesignType()
obj_dtype = my_moop.getObjectiveType()
sim_dtype = my_moop.getSimulationType()
# Display the dtypes as strings
print("Design variable type: " + str(des_dtype))
print("Simulation output type: " + str(sim_dtype))
print("Objective type: " + str(obj_dtype))
The above code produces the following output.
Design variable type: [('MyDes', '<f8'), ('MyCat', '<i4')]
Simulation output type: [('MySim', '<f8', (2,))]
Objective type: [('MyObj', '<f8')]
Unnamed Output Types
The following unnamed_var_ex.py code demonstrates ParMOO’s output datatype when the
MOOP object is defined using unnamed variables.
import numpy as np
from parmoo import MOOP
from parmoo.searches import LatinHypercube
from parmoo.surrogates import GaussRBF
from parmoo.acquisitions import UniformWeights
from parmoo.optimizers import LocalGPS
# Fix the random seed for reproducibility
np.random.seed(0)
my_moop = MOOP(LocalGPS)
# Define a simulation to use below
def sim_func(x):
return np.array([(x[0]) ** 2, (x[0] - 1.0) ** 2])
# Add a design variable, simulation, objective, and constraint, w/o name key
my_moop.addDesign({'des_type': "continuous",
'lb': 0.0, 'ub': 1.0})
my_moop.addSimulation({'m': 2,
'sim_func': sim_func,
'search': LatinHypercube,
'surrogate': GaussRBF,
'hyperparams': {'search_budget': 20}})
my_moop.addObjective({'obj_func': lambda x, s: sum(s)})
my_moop.addConstraint({'constraint': lambda x, s: 0.1 - x[0]})
# Extract numpy dtypes for all of this MOOP's inputs/outputs
des_dtype = my_moop.getDesignType()
sim_dtype = my_moop.getSimulationType()
obj_dtype = my_moop.getObjectiveType()
const_dtype = my_moop.getConstraintType()
# Display the dtypes as strings
print("Design variable type: " + str(des_dtype))
print("Simulation output type: " + str(sim_dtype))
print("Objective type: " + str(obj_dtype))
print("Constraint type: " + str(const_dtype))
print()
# Add one acquisition and solve with 0 iterations to initialize databases
my_moop.addAcquisition({'acquisition': UniformWeights})
my_moop.solve(0)
# Extract final objective and simulation databases
obj_db = my_moop.getObjectiveData()
sim_db = my_moop.getSimulationData()
# Print the objective database dtypes
print("Objective database keys: " + str([key for key in obj_db.keys()]))
for key in obj_db.keys():
print("\t'" + key + "'" + " dtype: " + str(obj_db[key].dtype))
print("\t'" + key + "'" + " shape: " + str(obj_db[key].shape))
print()
# Print the simulation database dtypes
print("Simulation database type: " + str(type(sim_db)))
print("Simulation database length: " + str(len(sim_db)))
for i, dbi in enumerate(sim_db):
print("\tsim_db[" + str(i) + "] database keys: " +
str([key for key in dbi.keys()]))
for key in dbi.keys():
print("\t\t'" + key + "'" + " dtype: " + str(dbi[key].dtype))
print("\t\t'" + key + "'" + " shape: " + str(dbi[key].shape))
The above code produces the following output.
Design variable type: ('<f8', (1,))
Simulation output type: ('<f8', (2,))
Objective type: ('<f8', (1,))
Constraint type: ('<f8', (1,))
Objective database keys: ['x_vals', 'f_vals', 'c_vals']
'x_vals' dtype: float64
'x_vals' shape: (20, 1)
'f_vals' dtype: float64
'f_vals' shape: (20, 1)
'c_vals' dtype: float64
'c_vals' shape: (20, 1)
Simulation database type: <class 'list'>
Simulation database length: 1
sim_db[0] database keys: ['x_vals', 's_vals']
'x_vals' dtype: float64
'x_vals' shape: (20, 1)
's_vals' dtype: float64
's_vals' shape: (20, 2)
Adding Precomputed Simulation Values
The following precomputed_data.py code demonstrates how to add a precomputed simulation output to ParMOO’s internal simulation databases.
import numpy as np
from parmoo import MOOP
from parmoo.searches import LatinHypercube
from parmoo.surrogates import GaussRBF
from parmoo.acquisitions import UniformWeights
from parmoo.optimizers import LocalGPS
my_moop = MOOP(LocalGPS)
my_moop.addDesign({'name': "x1",
'des_type': "continuous",
'lb': 0.0, 'ub': 1.0})
my_moop.addDesign({'name': "x2", 'des_type': "categorical",
'levels': 3})
def sim_func(x):
if x["x2"] == 0:
return np.array([(x["x1"] - 0.2) ** 2, (x["x1"] - 0.8) ** 2])
else:
return np.array([99.9, 99.9])
my_moop.addSimulation({'name': "MySim",
'm': 2,
'sim_func': sim_func,
'search': LatinHypercube,
'surrogate': GaussRBF,
'hyperparams': {'search_budget': 20}})
my_moop.addObjective({'name': "f1", 'obj_func': lambda x, s: s["MySim"][0]})
my_moop.addObjective({'name': "f2", 'obj_func': lambda x, s: s["MySim"][1]})
my_moop.addAcquisition({'acquisition': UniformWeights})
# Precompute one simulation value for demo
des_val = np.zeros(1, dtype=[("x1", float), ("x2", int)])[0]
sim_val = sim_func(des_val)
# Add the precomputed simulation value from above
my_moop.update_sim_db(des_val, sim_val, "MySim")
# Get and display initial database
sim_db = my_moop.getSimulationData()
print(sim_db)
The above code produces the following output.
{'MySim': array([(0., 0, [0.04, 0.64])],
dtype=[('x1', '<f8'), ('x2', '<i4'), ('out', '<f8', (2,))])}
Logging and Checkpointing
When solving a large or expensive problem, logging and checkpointing are recommended. The following code snippet shows how to solve the same problem from the Quickstart example, but with logging and checkpointing turned on. Then, another MOOP is created by reloading from the checkpoint file, and run for 1 extra iteration, in order to demonstrate how to reload from a saved checkpoint file.
Note how the lambda functions from Quickstart example have been explicitly defined below. This is because ParMOO reloads functions by name, which means that it does not support checkpointing for lambda functions and any named function must be defined with the same name in the global scope during reloading.
import numpy as np
from parmoo import MOOP
from parmoo.searches import LatinHypercube
from parmoo.surrogates import GaussRBF
from parmoo.acquisitions import UniformWeights
from parmoo.optimizers import LocalGPS
import logging
# Fix the random seed for reproducibility
np.random.seed(0)
# Create a new MOOP
my_moop = MOOP(LocalGPS)
# Add 1 continuous and 1 categorical design variable
my_moop.addDesign({'name': "x1",
'des_type': "continuous",
'lb': 0.0, 'ub': 1.0})
my_moop.addDesign({'name': "x2", 'des_type': "categorical",
'levels': 3})
# Create a simulation function
def sim_func(x):
if x["x2"] == 0:
return np.array([(x["x1"] - 0.2) ** 2, (x["x1"] - 0.8) ** 2])
else:
return np.array([99.9, 99.9])
# Add the simulation function to the MOOP
my_moop.addSimulation({'name': "MySim",
'm': 2,
'sim_func': sim_func,
'search': LatinHypercube,
'surrogate': GaussRBF,
'hyperparams': {'search_budget': 20}})
# Define the 2 objectives as named Python functions
def obj1(x, s): return s["MySim"][0]
def obj2(x, s): return s["MySim"][1]
# Define the constraint as a function
def const(x, s): return 0.1 - x["x1"]
# Add 2 objectives
my_moop.addObjective({'name': "f1", 'obj_func': obj1})
my_moop.addObjective({'name': "f2", 'obj_func': obj2})
# Add 1 constraint
my_moop.addConstraint({'name': "c1", 'constraint': const})
# Add 3 acquisition functions (generates batches of size 3)
for i in range(3):
my_moop.addAcquisition({'acquisition': UniformWeights,
'hyperparams': {}})
# Turn on logging with timestamps
logging.basicConfig(level=logging.INFO,
format='%(asctime)s %(levelname)-8s %(message)s',
datefmt='%Y-%m-%d %H:%M:%S')
# Use checkpointing without saving a separate data file (in "parmoo.moop" file)
my_moop.setCheckpoint(True, checkpoint_data=False, filename="parmoo")
# Solve the problem with 4 iterations
my_moop.solve(4)
# Create a new MOOP object and reload the MOOP from parmoo.moop file
new_moop = MOOP(LocalGPS)
new_moop.load("parmoo")
# Do another iteration
new_moop.solve(5)
# Display the solution
results = new_moop.getPF()
print(results, "\n dtype=" + str(results.dtype))
The above checkpointing.py code produces the following output.
[(0.79013666, 0, 3.48261275e-01, 9.72855201e-05, -0.69013666)
(0.7895263 , 0, 3.47541259e-01, 1.09698380e-04, -0.6895263 )
(0.73488156, 0, 2.86098283e-01, 4.24041125e-03, -0.63488156)
(0.70656124, 0, 2.56604287e-01, 8.73080244e-03, -0.60656124)
(0.68409101, 0, 2.34344111e-01, 1.34348928e-02, -0.58409101)
(0.67237225, 0, 2.23135544e-01, 1.62888423e-02, -0.57237225)
(0.5784217 , 0, 1.43202981e-01, 4.90969442e-02, -0.4784217 )
(0.53031761, 0, 1.09109723e-01, 7.27285920e-02, -0.43031761)
(0.51322778, 0, 9.81116425e-02, 8.22383058e-02, -0.41322778)
(0.49723345, 0, 8.83477213e-02, 9.16675863e-02, -0.39723345)
(0.44987019, 0, 6.24351129e-02, 1.22590882e-01, -0.34987019)
(0.40591372, 0, 4.24004606e-02, 1.55303995e-01, -0.30591372)
(0.39028872, 0, 3.62097978e-02, 1.67863331e-01, -0.29028872)
(0.38995793, 0, 3.60840145e-02, 1.68134501e-01, -0.28995793)
(0.38287784, 0, 3.34443041e-02, 1.73990897e-01, -0.28287784)
(0.25435646, 0, 2.95462529e-03, 2.97726867e-01, -0.15435646)
(0.20137796, 0, 1.89878434e-06, 3.58348342e-01, -0.10137796)]
dtype=[('x1', '<f8'), ('x2', '<i4'), ('f1', '<f8'), ('f2', '<f8'), ('c1', '<f8')]
Solving a MOOP with Derivative-Based Solvers
This example shows how to minimize two conflicting quadratic functions
of three variables (named x1, x2, and x3),
under the constraint that an additional categorical variable (x4)
must be fixed in class 0 \(^*\),
using the derivative-based solver LBFGSB.
\(^*\) No, this constraint does not really affect the solution;
it is just here to demonstrate how constraints/categorical variables
are handled by derivative-based solvers.
ParMOO does not use derivative information associated with
any categorical variables, so the derivative w.r.t. x4 can be set
to any value, without affecting the outcome.
import numpy as np
from parmoo import MOOP
from parmoo.acquisitions import RandomConstraint, FixedWeights
from parmoo.searches import LatinHypercube
from parmoo.surrogates import GaussRBF
from parmoo.optimizers import LBFGSB
# Fix the random seed for reproducibility
np.random.seed(0)
# Create a new MOOP with a derivative-based solver
my_moop = MOOP(LBFGSB, hyperparams={})
# Add 3 continuous variables named x1, x2, x3
for i in range(3):
my_moop.addDesign({'name': "x" + str(i+1),
'des_type': "continuous",
'lb': 0.0,
'ub': 1.0,
'des_tol': 1.0e-8})
# Add one categorical variable named x4
my_moop.addDesign({'name': "x4",
'des_type': "categorical",
'levels': 3})
def quad_sim(x):
""" A quadratic simulation function with 2 outputs.
Returns:
np.ndarray: simulation value (S) with 2 outputs
* S_1(x) = <x, x>
* S_2(x) = <x-1, x-1>
"""
return np.array([x["x1"] ** 2 + x["x2"] ** 2 + x["x3"] ** 2,
(x["x1"] - 1.0) ** 2 + (x["x2"] - 1.0) ** 2 +
(x["x3"] - 1.0) ** 2])
# Add the quadratic simulation to the problem
# Use a 10 point LH search for ex design and a Gaussian RBF surrogate model
my_moop.addSimulation({'name': "f_conv",
'm': 2,
'sim_func': quad_sim,
'search': LatinHypercube,
'surrogate': GaussRBF,
'hyperparams': {'search_budget': 10}})
def obj_f1(x, sim, der=0):
""" Minimize the first output from 'f_conv' """
if der == 0:
return sim['f_conv'][0]
elif der == 1:
return np.zeros(1, dtype=x.dtype)[0]
elif der == 2:
result = np.zeros(1, dtype=sim.dtype)[0]
result['f_conv'][0] = 1.0
return result
def obj_f2(x, sim, der=0):
""" Minimize the second output from 'f_conv' """
if der == 0:
return sim['f_conv'][1]
elif der == 1:
return np.zeros(1, dtype=x.dtype)[0]
elif der == 2:
result = np.zeros(1, dtype=sim.dtype)[0]
result['f_conv'][1] = 1.0
return result
# Minimize each of the 2 outputs from the quadratic simulation
my_moop.addObjective({'name': "f1",
'obj_func': obj_f1})
my_moop.addObjective({'name': "f2",
'obj_func': obj_f2})
def const_x4(x, sim, der=0):
""" Constrain x["x4"] = 0 """
if der == 0:
return 0.0 if (x["x4"] == 0) else 1.0
elif der == 1:
# No derivatives for categorical design var, just return all zeros
return np.zeros(1, dtype=x.dtype)[0]
elif der == 2:
return np.zeros(1, dtype=sim.dtype)[0]
# Add the single constraint to the problem
my_moop.addConstraint({'name': "c_x4",
'constraint': const_x4})
# Add 2 different acquisition functions to the problem
my_moop.addAcquisition({'acquisition': RandomConstraint})
my_moop.addAcquisition({'acquisition': FixedWeights,
# Fixed weight with equal weight on both objectives
'hyperparams': {'weights': np.array([0.5, 0.5])}})
# Turn on checkpointing -- creates the files parmoo.moop and parmoo.surrogate.1
my_moop.setCheckpoint(True, checkpoint_data=False, filename="parmoo")
# Turn on logging
import logging
logging.basicConfig(level=logging.INFO,
format='%(asctime)s %(levelname)-8s %(message)s',
datefmt='%Y-%m-%d %H:%M:%S')
# Solve the problem
my_moop.solve(5)
# Get and print full simulation database
sim_db = my_moop.getSimulationData(format="pandas")
print("Simulation data:")
for key in sim_db.keys():
print(f"\t{key}:")
print(sim_db[key])
# Get and print results
soln = my_moop.getPF(format="pandas")
print("\n\n")
print("Solution points:")
print(soln)
The above advanced_ex.py code produces the following output.
Simulation data:
f_conv:
x1 x2 x3 x4 out_0 out_1
0 0.164589 0.071519 0.256804 0 0.098153 2.112328
1 0.279173 0.761210 0.446148 2 0.856425 0.883365
2 0.054881 0.143759 0.645615 1 0.440497 1.751987
3 0.761764 0.514335 0.869763 1 1.601312 0.309588
4 0.563992 0.252889 0.383262 1 0.528930 1.128643
5 0.626456 0.302022 0.761693 2 1.063841 0.683499
6 0.835951 0.921038 0.912893 1 2.380498 0.040735
7 0.308713 0.677423 0.189177 0 0.589994 1.239367
8 0.497862 0.843703 0.060276 0 0.963335 1.159652
9 0.967064 0.479916 0.594467 0 1.518922 0.436029
10 0.198616 0.279706 0.253141 0 0.181764 1.718838
11 0.104311 0.713583 0.360565 0 0.650089 1.293170
12 0.990617 0.496152 0.621227 2 1.613413 0.397420
13 0.322478 0.475506 0.173748 1 0.360287 1.416822
14 0.143465 0.160993 0.217433 0 0.093778 2.049997
15 0.317939 0.476173 0.171654 0 0.357290 1.425761
16 0.711117 0.366636 0.640507 0 1.050358 0.613839
17 0.373639 0.410939 0.465101 0 0.524796 1.025438
18 0.193532 0.456200 0.200312 0 0.285698 1.585610
19 0.448389 0.368552 0.567048 0 0.658426 0.890449
Solution points:
x1 x2 x3 x4 f1 f2 c_x4
0 0.967064 0.479916 0.594467 0 1.518922 0.436029 0.0
1 0.711117 0.366636 0.640507 0 1.050358 0.613839 0.0
2 0.448389 0.368552 0.567048 0 0.658426 0.890449 0.0
3 0.373639 0.410939 0.465101 0 0.524796 1.025438 0.0
4 0.317939 0.476173 0.171654 0 0.357290 1.425761 0.0
5 0.193532 0.456200 0.200312 0 0.285698 1.585610 0.0
6 0.198616 0.279706 0.253141 0 0.181764 1.718838 0.0
7 0.143465 0.160993 0.217433 0 0.093778 2.049997 0.0
Note how in the full simulation database several of the design points
violate the constraint (x4 != 0).
But in the solution, the constraint is always satisfied.