18  Simulating an orchestrated platform

This chapter ties together several of the previous lessons to demonstrate how the process logic of an orchestrated platform can be simulated with PyLabRobot. This platform is meant to be a realistic example of a workcell in a highly automated lab handling stochastic inputs. This use case exemplifies many realistic demands that will be placed on an automated workcell - optimizing pipetting operations for time efficiency, tracking reagent use over time, recording all actions and data in an organized way.

Modeling how a platform responds to realistic data is invaluable for process design that can be extended to many different protocols. We encourage the reader to consider how this methodology could be applied to designing other processes and platforms.

Long waits and physical measurements are two reasons laboratory code is difficult to test quickly. This chapter substitutes simulated time and stochastic measurements, then combines them with the run directory of chapter 13, the SQLite tables of chapter 14, and the wrappers of chapter 15 into a single simulated unit operation: a normalization cell fed by a plate stacker.

NoteSetup

import asyncio, json, logging, random, sqlite3, time
from datetime import datetime
from pathlib import Path

from pylabrobot.liquid_handling import LiquidHandler
from pylabrobot.liquid_handling.backends import LiquidHandlerChatterboxBackend
from pylabrobot.resources import (
    STARLetDeck, PLT_CAR_L5AC_A00, cor_96_wellplate_360uL_Fb,
    nest_1_troughplate_195000uL_Vb, TIP_CAR_288_C00,
    opentrons_96_filtertiprack_200ul, ResourceStack, Coordinate,
    set_volume_tracking, set_tip_tracking,
)

---

18.1 Simulate time and data together

Simulated time is a highly useful concept for testing methods whose success is contingent on timing and throughput rate.

A clock that can run against either simulated or real time:

class Clock:
    def __init__(self, simulation: bool):
        self.simulation = simulation
        self.time = 0.0

    def now(self):
        if self.simulation:
            return self.time
        return time.monotonic()

    async def wait(self, seconds):
        if self.simulation:
            self.time += seconds
        else:
            await asyncio.sleep(seconds)

    async def wait_until(self, when):
        await self.wait(max(0.0, when - self.now()))

You can use this class to get simulated time from a simulated method, and real time during a live method, with just a single Boolean variable passed to the simulating argument.

The same idea applies to measurements. You can simulate credible distributions of random variables, like DNA concentrations in samples, and use these to test the robustness of your process against random data.

rng = random.Random(42)

def quantify(true_conc, cv=0.03):
    return true_conc * (
        1 + rng.gauss(0, cv)
    )

---

18.2 Unit operations

Normalizing concentrations in samples is one of the most widespread operations in lab automation. Here we will design a workcell that can handle an asynchronously populated stack of 96-well plates with samples that we must normalize. Arrival times are random, numbers of samples are random, concentrations are random.

arrivals ──► input stacker ──► liquid handler ──► output stacker
             (plate hotel)      (dilute to target)

Let’s identify the main unit operations we need and break them down into functions. This is one of the most valuable architectural parts of scripting a protocol. Not only will this make the final script more readable, like a natural language protocol, but identifying unit operations will make it easier to reason about the logical flow of your process and will make each operation individually testable with defined inputs and outputs.

We don’t cover testing strategy in this cookbook, since it would take a very large amount of content to do justice to the topic, and it is highly contingent on physical implementation details outside of our scope. Suffice to say that a robust platform is built from rigorously tested unit operations that are composed into further unit operations and tested at each level of composition. The root causes of failures, as they inevitably occur, are identified and addressed through tests on unit operations. In a sense, the test is the fundamental unit and measure of certainty in an automated platform, just as an experiment is the measure of certainty for a scientific hypothesis.

18.2.1 Build the cell

The STAR’s iSWAP moves the plates, so transport and pipetting contend for one instrument: the arm cannot fetch the next plate while the head is dispensing. Everything else in the cell — carrier, trough, tips, both stackers — is passive labware with no backend behind it.

lh = LiquidHandler(
    backend=LiquidHandlerChatterboxBackend(),
    deck=STARLetDeck(),
)
await lh.setup()

carrier = PLT_CAR_L5AC_A00(name="carrier")
lh.deck.assign_child_resource(carrier, rails=10)
carrier[0] = buffer = nest_1_troughplate_195000uL_Vb(
    name="buffer"
)

tip_carrier = TIP_CAR_288_C00(name="tip_carrier")
lh.deck.assign_child_resource(tip_carrier, rails=2)
tip_carrier[0] = rack = opentrons_96_filtertiprack_200ul(
    name="rack"
)

The two stackers here are ResourceStacks, which are just convenience objects for stacks of labware.

inp = ResourceStack(name="input_stacker", direction="z")
out = ResourceStack(name="output_stacker", direction="z")

lh.deck.assign_child_resource(
    inp, location=Coordinate(600, 300, 100)
)
lh.deck.assign_child_resource(
    out, location=Coordinate(600, 100, 100)
)

set_volume_tracking(True)
set_tip_tracking(True)
buffer["A1"][0].set_volume(195_000)

That is one rack of 96 filtered 200 µL tips and a single-well trough holding 195 mL.

18.2.2 Dilution calculation

Each well holds V_SAMPLE µL at a known concentration. Diluting to TARGET means adding buffer until the volume has grown by the ratio of the two concentrations.

TARGET, V_SAMPLE, WELL_MAX = 10.0, 10.0, 360.0

def plan_well(conc):
    if conc <= TARGET:
        return 0.0, "under"

    buffer_vol = V_SAMPLE * (conc / TARGET - 1)

    if V_SAMPLE + buffer_vol > WELL_MAX:
        return 0.0, "overflow"

    return round(buffer_vol, 1), "ok"

18.2.3 Record the run

The run directory and the process log are chapter 13, as one function:

def new_run(protocol, machines):
    run_id = datetime.now().strftime("%Y%m%d_%H%M%S")
    run_dir = Path("experiments") / protocol / f"{run_id}_{protocol}"

    for subdir in ["inputs", "raw", "derived", "logs", "state"]:
        (run_dir / subdir).mkdir(parents=True, exist_ok=True)

    manifest = {
        "run_id": run_id,
        "started": datetime.now().isoformat(),
        "protocol": protocol,
        "plates": [],
        "machines": machines,
    }
    (run_dir / "manifest.json").write_text(
        json.dumps(manifest, indent=2)
    )

    log = logging.getLogger(f"run.{run_id}")
    log.setLevel(logging.INFO)
    log.propagate = False
    log.addHandler(
        logging.FileHandler(run_dir / "logs" / "process.log")
    )

    return run_dir, manifest, log

1

Filled in at the end of the run, once the cell knows which plates it actually saw. The sample flow is stochastic, so no manifest written at the start can name them.

2

The process log is this run’s file, not the root logger’s console. Without this, every line would also land wherever the surrounding application happens to log.

The tables are chapter 14 — one row per plate, one row per well:

def open_tables(path):
    db = sqlite3.connect(path)

    db.execute("""
        CREATE TABLE plates (
            plate     TEXT PRIMARY KEY,
            arrived   REAL,
            concs     TEXT,
            started   REAL,
            finished  REAL
        )
    """)

    db.execute("""
        CREATE TABLE wells (
            plate     TEXT,
            well      INTEGER,
            conc      REAL,
            buffer    REAL,
            flag      TEXT,
            achieved  REAL
        )
    """)

    return db

Which makes starting a run two lines:

run_dir, manifest, log = new_run(
    "normalization", machines=["liquid_handler"]
)
db = open_tables(run_dir / "run.sqlite")

18.2.4 Model the arrivals

Arrivals are the sample flow into this cell, reduced to the only two things the cell needs from it: when a plate shows up, and what is in it. The sampling is a generator — suspended between arrivals, holding t and n as ordinary local variables, drawing a plate only when asked:

sim_rng = random.Random(42)
MEAN_GAP, SERVICE, WINDOW, POLL = 120, 300, 900, 30

def arrivals(rng, start, until):
    t, n = start, 0

    while True:
        t += rng.expovariate(1 / MEAN_GAP)

        if t > until:
            return

        n += 1
        yield t, f"plate_{n}", [
            max(1.0, rng.gauss(60, 25))
            for _ in range(rng.randint(3, 8))
        ]

1

Plates arrive faster than the cell can process them — a mean gap of 120 s against a 300 s service time. That is the interesting case; with a gap longer than the service time the stacker never holds more than one plate and nothing below is visible. WINDOW is how long samples keep arriving; POLL is how often the cell looks.

2

Arrivals start from the clock’s own zero, so the same code is valid whether now() begins at virtual zero or at the operating system’s current monotonic value.

3

The gap is drawn before the cutoff is checked, so the window closes on a real draw rather than on a plate count decided in advance.

4

A plate off an extraction step is rarely full. Varying the count is what keeps the channel arithmetic below honest: nothing may assume eight.

An exponential gap means the process has no memory: how long the cell has already waited says nothing about how much longer it will. That is what makes the stacker depth, and the wait times below, worth simulating rather than reasoning about.

18.2.5 Simulated and real sample flows

The orchestration loop must not be written against the generator. A simulated sample flow can be asked when the next plate will arrive; a real one cannot. Build the loop on the two questions reality can answer — what has arrived by now, and has the flow stopped — and both sides fit:

class SimulatedSampleFlow:
    def __init__(self, rng, start, until):
        self.stream = arrivals(rng, start, until)
        self.pending = next(self.stream, None)

    def due(self, now):
        ready = []

        while self.pending is not None and self.pending[0] <= now:
            ready.append(self.pending)
            self.pending = next(self.stream, None)

        return ready

    def closed(self, now):
        return self.pending is None

1

The source may know the future — that is what makes it a simulation. The loop may not, which is what keeps the loop runnable on hardware.

2

Every plate that has appeared since the last call: none, one, or several if the cell stayed busy through more than one gap.

3

No more samples are coming. The loop drains what is left and exits.

The real implementation answers the same two questions by looking rather than by drawing. Write this one against your own barcode reader and sample database; the loop does not change:

#| eval: false
class StackerSampleFlow:
    def __init__(self, reader, samples, until):
        self.reader, self.samples = reader, samples
        self.until = until
        self.seen = set()

    def due(self, now):
        scanned = self.reader.scan()
        fresh = [p for p in scanned if p not in self.seen]
        self.seen.update(fresh)

        return [
            (now, p, self.samples.concentrations(p))
            for p in fresh
        ]

    def closed(self, now):
        return now > self.until and not self.reader.scan()

1

now is when the cell noticed the plate, not when it was set down — the only arrival time a real cell ever has. In simulation the two are identical, which is one of the things simulation quietly makes easier than reality.

18.2.6 Load an arrived plate

Each arrived plate is set on top of the input stacker by whatever feeds this cell. The robot is not involved, and the row is written at the moment the plate is real — the database records what happened, never what is scheduled to. This function is shared by both sample flows: on hardware the plate is already sitting there, and assign_child_resource is how the resource tree is told about a plate that exists:

def load_plate(arrived, plate_id, concs):
    plate = cor_96_wellplate_360uL_Fb(name=plate_id)

    for well in plate[: len(concs)]:
        well.set_volume(V_SAMPLE)

    inp.assign_child_resource(plate)
    db.execute(
        "INSERT INTO plates(plate, arrived, concs) VALUES (?,?,?)",
        (plate_id, arrived, json.dumps(concs)),
    )
    db.commit()

1

plate[:n] is the integer-slice form of the selection grammar — the first n wells in fill order, already unwrapped, so the body reads well.set_volume(...) rather than indexing back into the plate. Only the occupied wells get a volume; the rest are empty, and the volume tracker will say so if the protocol reaches for one.

18.2.7 Normalize a plate

handle_errors is the wrapper from chapter 15. It is retyped here rather than imported, so this chapter runs on its own:

import functools

def handle_errors(fn):
    @functools.wraps(fn)
    async def wrapper(*args, **kwargs):
        try:
            return await fn(*args, **kwargs)
        except Exception:
            log.exception("%s failed", fn.__name__)
            raise
    return wrapper

It writes to the run’s own process log, so a failure is recorded in the run directory next to the data it failed to produce.

Now one plate through the cell. Two things vary per plate — how many samples it carries, and which of those need buffer at all — and both are handled the same way, by choosing channels rather than by padding volumes:

@handle_errors
async def normalize(plate, concs, tip_column):
    rows = [plan_well(c) for c in concs]
    need = [i for i, (vol, _) in enumerate(rows) if vol > 0]
    vols = [rows[i][0] for i in need]

    if not need:
        return rows

    await lh.pick_up_tips(
        [tip_column[i] for i in need], use_channels=need
    )
    await lh.aspirate(
        [buffer["A1"][0]] * len(need),
        vols=vols,
        use_channels=need,
    )
    await lh.dispense(
        [plate[i][0] for i in need],
        vols=vols,
        use_channels=need,
    )
    await lh.discard_tips(use_channels=need)

    return rows

1

The plate’s sample count comes from the worklist, not from the labware — concs is however many wells were filled upstream, and everything below is sized from it.

2

use_channels must match the length of vols; a well that needs nothing is simply not a channel in this operation. Between a partial plate and the wells already at target, the channels in use are a subset of a subset. The same indices select the tips: channel i takes spot i of the plate’s tip column.

3

A plate whose wells are all at or below target needs no buffer at all. Picking up zero tips is not a no-op worth asking the robot to perform.

That is one aspirate and one dispense for the whole plate, not a loop of single transfers. The backend sees every active channel move together, each carrying its own volume: use_channels=need says which channels act, no channel is pinned to a rack row, and the trough offsets are not asked for — PLR spreads the active channels across the reservoir on its own.

18.2.8 Read and write the run’s tables

Everything the run reads or writes is behind a name, so the loop itself is only the protocol:

def worklist(plate):
    (concs,) = db.execute(
        "SELECT concs FROM plates WHERE plate=?", (plate.name,)
    ).fetchone()
    return json.loads(concs)

def record_plate(plate, samples, rows, started, finished):
    for i, ((vol, flag), conc) in enumerate(zip(rows, samples)):
        db.execute(
            "INSERT INTO wells VALUES (?,?,?,?,?,?)",
            (
                plate.name, i, conc, vol, flag,
                quantify(conc * V_SAMPLE / (V_SAMPLE + vol)),
            ),
        )

    (arrived,) = db.execute(
        "SELECT arrived FROM plates WHERE plate=?", (plate.name,)
    ).fetchone()
    db.execute(
        "UPDATE plates SET started=?, finished=? WHERE plate=?",
        (started, finished, plate.name),
    )
    db.commit()
    log.info("%s waited %.0f s", plate.name, started - arrived)

def close_run():
    manifest["plates"] = [
        plate_id for (plate_id,) in db.execute(
            "SELECT plate FROM plates ORDER BY arrived"
        )
    ]
    manifest["finished"] = datetime.now().isoformat()
    (run_dir / "manifest.json").write_text(
        json.dumps(manifest, indent=2)
    )

1

The achieved concentration is what the dilution actually produced, seen through a quantification carrying the noise from the first recipe.

---

18.3 Run the cell

18.3.1 Pick simulation or hardware

One flag picks both the clock and the sample flow, and nothing below it knows which it got:

SIMULATION = True

clock = Clock(simulation=SIMULATION)
sample_flow = (
    SimulatedSampleFlow(sim_rng, start=clock.now(), until=WINDOW)
    if SIMULATION
    else StackerSampleFlow(reader, samples, until=WINDOW)
)

1

The whole switch. Simulated, this chapter renders in about a second; real, the same loop paces itself against elapsed time and a barcode reader.

2

Only the branch that is taken is evaluated, so this cell runs with StackerSampleFlow undefined — as it is here, where there is no reader to hand it.

18.3.2 Run the loop

The loop has no scheduler in it. It admits whatever arrived, takes the top plate if there is one, and otherwise waits one poll interval. The hardware chooses which plate runs, and the loop only records the consequence:

tip_columns = (rack.column(c) for c in range(rack.num_items_x))

while True:
    for arrival in sample_flow.due(clock.now()):
        load_plate(*arrival)

    if len(inp.children) == 0:
        if sample_flow.closed(clock.now()):
            break

        await clock.wait(POLL)
        continue

    plate = inp.get_top_item()
    samples = worklist(plate)
    started = clock.now()

    await lh.move_plate(plate, carrier[1])
    rows = await normalize(plate, samples, next(tip_columns))
    await lh.move_plate(plate, out)
    await clock.wait(SERVICE)

    record_plate(plate, samples, rows, started, clock.now())

close_run()

1

One tip column per plate — eight spots, one per channel, out of the rack’s twelve columns. A plate off an extraction step carries at most eight samples, so a column covers any of them, and the plate takes only the spots its use_channels names. When the rack is spent, next() raises StopIteration, which is the true statement: this cell has run out of tips.

2

The only question the loop asks about arrivals, and one a real cell can answer.

3

Empty stacker and no more samples coming: everything that arrived has been processed.

4

Idle only because the stacker is empty and nothing has arrived. Simulated, the poll costs nothing and the clock jumps POLL seconds; real, it is a POLL-second sleep. The loop cannot tell.

5

get_top_item() is the only plate the gripper can reach.

6

A plate’s tips are its own column, so nothing is shared between plates and the mapping from spot to channel is the same index in both.

for plate_id, arrived, started in db.execute(
    "SELECT plate, arrived, started FROM plates ORDER BY arrived"
):
    print(
        f"{plate_id}  arrived {arrived:6.0f}  "
        f"started {started:6.0f}  waited {started - arrived:6.0f}"
    )

plate_1  arrived    122  started    150  waited     28
plate_2  arrived    284  started   1950  waited   1666
plate_3  arrived    311  started   1650  waited   1339
plate_4  arrived    350  started    450  waited    100
plate_5  arrived    576  started    750  waited    174
plate_6  arrived    836  started   1350  waited    514
plate_7  arrived    877  started   1050  waited    173

18.3.3 Test process invariants

The exact numbers depend on the seeded arrivals; these properties do not.

Every plate was processed, and none before it arrived:

assert db.execute(
    "SELECT COUNT(*) FROM plates WHERE finished IS NULL"
).fetchone()[0] == 0

assert db.execute(
    "SELECT COUNT(*) FROM plates WHERE started < arrived"
).fetchone()[0] == 0
print("every plate processed, none before it arrived")

every plate processed, none before it arrived

Every sample that arrived has exactly one well record — no plate was half-processed, and no plate was sized against the labware instead of its worklist:

arrived_samples = sum(
    len(json.loads(concs))
    for (concs,) in db.execute("SELECT concs FROM plates")
)
recorded = db.execute("SELECT COUNT(*) FROM wells").fetchone()[0]

assert arrived_samples == recorded
print(f"{recorded} samples arrived, {recorded} recorded")

print("per plate:", [
    n for (n,) in db.execute(
        "SELECT COUNT(*) FROM wells GROUP BY plate ORDER BY plate"
    )
])

37 samples arrived, 37 recorded
per plate: [3, 7, 8, 4, 7, 5, 3]

Every well reached target or was flagged:

for conc, vol, flag, achieved in db.execute(
    "SELECT conc, buffer, flag, achieved FROM wells"
):
    if flag == "ok":
        assert abs(achieved - TARGET) / TARGET <= 0.10
print("every unflagged well within 10% of target")

every unflagged well within 10% of target

18.3.4 Exercise: how old does the oldest plate get?

The run recorded arrived and started for every plate but never asked the question that matters. Nothing here depends on the average wait; what matters is whether the oldest plate on the deck was still good when it was processed.

#| eval: false
db.execute("SELECT MAX(started - arrived) FROM plates").fetchone()[0]

That is one number from one run — a single draw. The simulation is seeded and costs no wall-clock time, which is exactly what makes it cheap to ask for the distribution instead: re-run across seeds, collect the maximum from each, and look at the spread. Two things to watch for:

• The mean wait will not move between LIFO and FIFO. The maximum will, by a lot. Whichever statistic is measured decides which conclusion is reached.
• The parameter governing both is the ratio of SERVICE to MEAN_GAP. Push it toward 1 and see which statistic degrades first.

If the tail is longer than the chemistry tolerates, the fix is not a better query. It is different equipment.

---

18.4 What to remember

• Put the simulation switch at the boundaries, not through the protocol: one flag picks the clock and the sample flow, and the orchestration loop between them is the same code either way.
• The loop may only ask questions reality can answer. “What has arrived by now?” has an answer on hardware; “when will the next plate arrive?” does not. One lookahead is all it takes to make a loop that only ever runs in simulation.
• Substitute time and measurements: a protocol that depends on a clock object and a seeded RNG can run instantly and reproduce a failing test.
• Seed the RNG (random.Random(42)) — stochastic, but deterministic.
• Keep the arithmetic that can be wrong in a plain function of numbers, separate from the robot, the database, and the clock.
• Test invariants — completion, ordering, tolerance to target, and physical position — rather than individual values. Set the tolerance wider than the simulated measurement noise.
• The resource tree is the check that a move was physically possible; chatterbox validates state, not physics.
• Sample the incoming flow inside the loop, one arrival at a time, against the same clock the cell runs on. A schedule generated in advance is a replay, not a simulation of arrivals.
• Polling is what makes the switch work on both sides: wait(POLL) is a free jump in simulation and a real sleep on hardware, and it costs the run only that a start time is rounded up to the next poll.
• Queue discipline is a property of the equipment. A ResourceStack is a stack, and serving from it is LCFS: same mean wait as FIFO, much heavier tail.