OpenEMS Optimizing Simulator
Introduction
Over at EOI we’ve been working on a time domain reflectometer (TDR) for precision agriculture for the past few months. The waveguide assembly has several competing design constraints. The following is an overview of a a surprisingly effective parameterization and optimization strategy for designing the waveguides structure. Using OpenEMS, OpenSCAD, SciPy, Docker, and a bit of Python we had an optimized form factor after a few days.
Here’s an example of the resulting waveform, showing a 214 ps rise time.
Here’s what we did.
System Overview
This optimization system includes:
- A parameterized OpenSCAD model for geometry generation.
- A Python simulation that uses Docker to run OpenEMS. Parameters are passed in using a CLI.
- A scoring function that evaluates TDR-like pulse responses.
- An optimizer from SciPy to minimize rise time and other signal artifacts.
- A STL generation step which generates the models from the OpenSCAD source.
- Simulation results stored in timestamped directories for later review.
flowchart LR
A[OpenSCAD Parameters] --> B[STL]
B --> C[Docker/OpenEMS Simulation]
C --> D[Pulse Response]
D --> E[Score]
E --> F[Optimizer]
F -->|Update parameters| ASimulation and optimization loop structure.
Running a script
docker run -v $(pwd):/home/appuser -t --rm --name openems openems python3 simulation_script.py --maxtime 10e-9
simulation/
├── simulation_script.py
├── caller_cli.py
├── objects/
| └── ...
├── results/
│ └── {RUN_NAME}_{VARIABLES}_{TIMESTAMP}/
│ ├── et
│ ├── ht
│ ├── port_it_1
│ ├── port_ut_1
│ ├── ...
│ ├── PEC_dump.vtp
│ ├── settings.csv
Calling OpenEMS from a script
The following executes the simulation script and passes in parameters as command line arguments. Use this to step through different parameters, in the next step we’ll use this to setup an optimizer.
def run_simulation(dname, script, kwargs):
args_ = []
kwargs.update({
"dir": str(dname),
"stl": str(dname / STL_FILE)})
for key, value in kwargs.items():
args_.append(f"--{key}")
args_.append(str(value))
args = ["docker", "run","-v", f"{base_dir_}:/home/appuser", "-t", "--rm", "--name", "dummy", "openems", "python3", str(script)]
args.extend(args_)
with subprocess.Popen(args) as process:
process.wait()
Adding an optimizer
To optimize we need a scalar score for each run.
Run Score
In our example this is the minimizing the rise time, overshoot, and pulse artifacts of a TDR return. I calculated these using the control library though I’ve also been developing a more tailored alternative here.
import control
import numpy as np
def response_score(x, y, *, bounds=None, weights=None, vdiff_min=0.1):
if not bounds:
bounds=(1800, 6000)
if not weights:
weights = {"1/RiseTime": 2000, "Overshoot": 500, "Undershoot": 500, "VDiff": 1000}
mask = (x > min(bounds)) & (x < max(bounds))
x_bound = x[mask]
y_bound = y[mask]
# Ensure we have data after bounding
if len(x_bound) == 0 or len(y_bound) == 0:
raise ValueError("No data in the specified bounds.")
step_info = control.step_info(sysdata=(y_bound, T=x_bound))
values = {
"1/RiseTime": 1/step_info["RiseTime"] if step_info["RiseTime"] > 0 else 0,
"Overshoot": step_info["Overshoot"],
"Undershoot": step_info["Undershoot"],
"VDiff": step_info["SettlingMax"]-step_info["SettlingMin"]}
score = sum([values[key]*weights.get(key, 0) for key in values])
valid = (np.max(y_bound) - np.min(y_bound)) > vdiff_min
return int(valid)*score, values, weights
For this project I added a few other options that also optimized on some of the parameters to make the waveguides easier to make. This is easily added so I won’t complicate the example with it.
Choosing an optimization strategy
This usually takes a few steps. For this project I did a rough estimate of the structure using two parallel cylinders. I used those values as the starting values and then ran the simulation for constant steps around those values. This results in 81 simulations (4 parameters with 3 variations: 3⁴ = 81). How long this takes depends on the simulation size but don’t wait around.
import pandas as pd
settings = pd.DataFrame([
{"name": "dy", "start": 12.7, "min": 12.7, "max": 30},
{"name": "top_angle", "start": 50, "min": 25, "max": 50},
{"name": "body_length", "start": 25, "min": 12.7, "max": 25},
{"name": "base_diameter", "start": 19, "min": 12, "max": 3*INCH/4},
])
# Number of steps for each parameter (can be customized)
num_steps = 3
# Generate value ranges for each parameter
param_names = settings["name"].tolist()
value_grids = [
np.linspace(row["min"], row["max"], num_steps)
for _, row in settings.iterrows()
]
parameter_settings = [
dict(zip(param_names, combo))
for combo in itertools.product(*value_grids)
]
for parameter_setting in parameter_settings:
optimize_geometry(args=None, settings=parameter_setting)
After some analysis and eyeballing I chose the values for the optimizer which were tweaked a few more times.
Scipy has a bunch of optimization algorithms, I got the best results using dual_annealing for global optimization.
from scipy.optimize import dual_annealing
x0 = np.array(settings["start"])
bounds = [(line["min"], line["max"]) for _, line in settings.iterrows()]
dual_annealing(optimize_geometry, bounds=bounds,
args=(settings,),
maxiter=1000,
initial_temp=5230.0,
restart_temp_ratio=2e-2,
visit=2.62,
accept=-5.0,
maxfun=10000000.0,
x0=x0)
While iterating I got good results using minimize with the Powell method.
from scipy.optimize import minimize
x0 = np.array(settings["start"])
bounds = [(line["min"], line["max"]) for _, line in settings.iterrows()]
res = minimize(
optimize_geometry, x0, tol=1e-5,
method="powell",
args=(settings), bounds=bounds,
options={"maxiter": 10000})
Let’s now define optimize_geometry, which ties everything together.
Optimization function (optimize_geometry)
This function needs to match the interface expected by scipy and return a scalar value. For reproducibility and debugging, each run stores both the simulation output and its input parameters in a timestamped directory. This lets us go back and analyze all the runs for debugging and for growing some intuition. We also store the weights and the measured values in results.csv.
import time
from pathlib import Path
SCRIPT=Path("simulation_script.py")
sub_dir_ = Path("results").absolute()
def optimize_geometry(args, settings):
default_settings = {
"top_angle": 49.5,
"side_angle": SIDE_ANGLE,
"dz": 1*INCH,
"dy": 0.93*INCH,
"base_diameter": 19,
"body_length": 25,
}
kwargs = default_settings
kwargs.update({pt: args[i] for i, pt in enumerate(settings["name"])})
dname = f"simulation_{int(time.time())}"
dname_local = sub_dir_ / dname
os.mkdir(dname_local)
make_model(dname_local, **kwargs)
df = pd.DataFrame([{"name": key, "value": value} for key, value in kwargs.items()])
df.to_csv(dname_local / "settings.csv")
dname_docker = Path("/home/appuser") / sub_dir_.stem / dname
run_simulation(dname_docker, script=SCRIPT, kwargs=kwargs)
xv1, yv1 = load_data(dname_local, port="port_ut_1")
score, values, weights = response_score(xv1, yv1)
df = pd.DataFrame({"values": values, "weights": weights})
df.to_csv(dname_local / "results.csv")
return -1*score
Note The score is made negative as we are minimizing and our weights are positive.
Make Model
Here I’m using OpenSCAD to generate an STL, which is then used in the simulation. The OpenSCAD model and the simulation connection logic are outside the scope of this post, but here’s the python for a generic set of arguments. You can also do this using JSON if preferred.
The STL geometry is generated with a parameterized OpenSCAD script that models the coupling structure.
def make_model(dname, source_file, stl_file, **kwargs):
output = dname / stl_file
arguments = " ".join([f"-D {key}={value}" for key, value in kwargs.items()])
os.system(f"openscad {str(source_file)} -o {output} {arguments}")
Simulation Script
To keep this post focused, I won’t dive deep into simulation setup details but I’d like to discuss the interface I’m using. I have OpenEMS in a Docker container (You can find the Docker build here: openEMS Docker container) and the tools I’m using are local so this is a little strange. I pass in the arguments through command line arguments using the click library.
@click.option("--dir", "dir_", default=None)
@click.option("-v", type=int, default=0)
@click.option("--zcut", type=float, default=18*inch, help="Make the section above this the surrounding material, chopping off the tines")
@click.option("--zmax", type=float, default=250)
@click.option("--fields", is_flag=True, default=False)
@click.option("--run", type=str, default="simulation_output")
@click.option("--maxtime", type=float, default=10e-9)
@click.option("--e_abs", type=float, default=3.2, help="Dielectic constant of epoxy")
@click.option("--e_dirt", type=float, default=7, help="Dielectic constant of dirt")
@click.option("--k_dirt", type=float, default=0, help="Conductivity of dirt")
@click.option("--stl", type=str, default=os.path.join(os.getcwd(),'./mechanical/pyramid-coupling-with-anvil.stl'))
@click.option("--no_csx", is_flag=True, default=False)
@click.command()
def main(dir_, v, zcut, zmax, fields, run, maxtime, e_abs, e_dirt, k_dirt, stl, no_csx):
...
if __name__ == "__main__":
# Change current path to script file folder
abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
main()
Summary
This post glossed over the simulation specifics, but it demonstrates how even a small amount of code can drive effective geometry optimization using Docker, OpenEMS, and Python.
This post outlined a lightweight but effective framework for geometry optimization using OpenEMS. By combining:
- Parametric geometry generation with OpenSCAD
- Dockerized simulation runs for clean reproducibility
- Python-based score functions for signal quality metrics
- And global and local optimization via SciPy
A more detailed walkthrough with a sharable version is in the works—stay tuned!