Heater Simulation using synApps epid and transform records#

TODO

Needs GUI screen shots, collected data, …

pushd /tmp/docker_ioc/iocgp/tmp/screens/ui/
caQtDM -macro "P=gp:,T=userTran1" yyTransform_full.ui &
caQtDM -macro "P=gp:,PID=epid1,TITLE=fb_epid" ./pid_control.ui &
popd

Simulate a temperature controller with heater and cooling features. Use it as a positioner.

A temperature sensor can be simulated with a single synApps swait record. The swait record has 16 fields for input values, either as constants or values from other EPICS PVs. Via a custom calculation expression, a simulated temperature sensor value is computed.

The simulation has these fields:

field	description
A	minimum “temperature” allowed
B	cooling rate parameter
C	heater power
D	output of PID loop
F	current “temperature”

A cooling simulation is necessary or the controller could never decrease the temperature.

This simulation is an example of building a custom ophyd Device.

Overview#

A PID loop (see below) has been paired with the simulated temperature using the synApps epid record, to update the operating power of the heater so that the temperature reaches a desired value.

Significant additional realism is added to the simulation by switching to the synApps transform record. The transform record also has 16 fields for values plus additional features:

value (same as swait record)
(optional) input PV link (same as swait record)
(optional) calculation expression
(optional) output PV link
(optional) descriptive text comment

Replacing swait with transform allows the addition of:

a random noise simulation
a tolerance value for evaluating if the setpoint and readback are in agreement
an at temperature feature

Schematic#

The PID loop decreases the difference (the following error) between the current temperature (readback) and a given setpoint by controlling the heater’s power:

power_fraction --> temperature --> PID --> output
  ^                                          |
  |                                          v
   <-----------------------------------------

Implementation#

The simulated temperature signal is computed from the various inputs of the transform record. The computation follows:

\(T_{sim} = T_{last} + \Delta T_{cool} + \Delta T_{heat} + \Delta T_{noise}\)

Certain fields of this transform record are linked to fields in an epid record. When the transform record is evaluated (processed, in EPICS terms), the epid is then processed. The epid record computes a new value (of power_fraction, the fraction of total heater power) for the next computation of the temperature signal.

The simulated temperature is re-computed (processed) periodically, as configured by the transform record’s .SCAN field. When transform is processed, it first pulls the power_fraction from epid (the output of the PID loop calculation) for the next simulated temperature.

Once transform processing is complete, the epid record is told to process itself by configuring the transform record’s forward link (.FLNK) field. The setpoint is pushed from transform to epid_ while epid pulls the current temperature (readback) from transform.

In addition to the simulated heater, the temperature calculation includes a simulation of cooling. This cooling is computed as a fraction of the last simulated temperature. Without any applied heating (when the heater power is OFF), the simulated will progress to the minimum temperature. (This simulation is a heater where the simulated temperature is constrained by T_min <= T <= T_max.)

PID Control Loops#

See these on-line references which explain PID control loops in detail.

ophyd interface#

This notebook will configure both transform and epid records using ophyd Python tools within the Bluesky framework. GUI screens from the caQtDM application will show the summary settings of each record.

record	EPICS PV	Python object	ophyd support
`transform`	`"gp:userTran1"`	`heater`	apstools.synApps.TransformRecord
`epid`	`"gp:epid1"`	`pid`	apstools.synApps.EpidRecord

Setup the `heater` and `pid`#

There’s a lot to configure here. Symbols have been defined as shortcuts to the transform record channels and comments have been added to explain.

The EpidRecord and TransformRecord classes connect with many PVs, far more than we need as an EPICS client to operate a temperature controller. We must access some of these fields to configure the records for our simulator. Local instances are created for each, within the setup function, to access all the settings for the EPICS configuration of the simulator. Once set, a simpler interface to the simulator PVs will be constructed.

Create a setup() which can be called when needed to reset the simulated heater.

[1]:

from apstools.synApps import EpidRecord
from apstools.synApps import TransformRecord

def heater_simulator_setup(transform_pv, epid_pv, starting=28.5, title="heater simulator"):
    """Configure the transform record as a heater with epid control."""

    transform = TransformRecord(transform_pv, name="transform")
    epid = EpidRecord(epid_pv, name="epid")
    transform.wait_for_connection()
    epid.wait_for_connection()

    transform.reset()  # clear all the transform record's configuration
    # TODO: epid.reset()  # no such thing available (now)

    transform.description.put(title)
    epid.description.put(f"{title} PID")
    transform.precision.put(3)  # note: real T only needs 1 sig fig

    # assign (local) Python symbols since channel names are long.
    # The calculation expressions (evaluated in the IOC) use channel letters.
    t_max = transform.channels.A
    t_min = transform.channels.B
    tolerance = transform.channels.C
    cooling = transform.channels.D
    power_fraction = transform.channels.E
    on_off = transform.channels.F
    noise = transform.channels.G
    t_last = transform.channels.H
    t_cooling = transform.channels.I
    t_heating = transform.channels.J
    t_noise = transform.channels.K
    smoothing = transform.channels.L
    setpoint = transform.channels.M
    readback = transform.channels.N
    difference = transform.channels.O
    at_temperature = transform.channels.P

    # simulated T will not go higher than this number
    # also used by heating & cooling
    t_max.comment.put("T max")
    t_max.current_value.put(500)

    # simulated T will not go lower than this number
    # also used by heating & cooling
    t_min.comment.put("T min")
    t_min.current_value.put(-10)

    # Acceptable range for difference between readback and setpoint
    tolerance.comment.put("tolerance")
    tolerance.current_value.put(1.0)  # same units as temperature

    # fraction to cool previous simulated temperature
    cooling.comment.put("cooling")
    cooling.current_value.put(0.05)  # 0.0 .. 1.0

    # PID will control this number from 0 (no power) to 1 (full power)
    power_fraction.comment.put("power fraction")
    power_fraction.current_value.put(0)  # 0.0 .. 1.0
    # assume EPICS IOC will take care of adding " NPP NMS" to the link
    power_fraction.input_pv.put(epid.output_value.pvname) # _from_ epid

    # User controls this ON/OFF switch.
    # No heating if the power is off.
    # This value will be passed _to_ the epid record.
    on_off.comment.put("heater ON")
    on_off.current_value.put(0)  # 0 or 1 (as float)
    on_off.output_pv.put(epid.feedback_on.pvname)  # _to_ epid

    # a random noise amplitude
    noise.comment.put("noise")
    noise.current_value.put(0.15)

    # Basis for the next computed temperature.
    # Previous computed temperature, smoothed by L.
    # Smoothing added to _reduce_ effects of simulated sensor noise.
    t_last.comment.put("T last")
    t_last.expression.put("H*L + N*(1-L)")  # 0.0 .. 1.0

    # temperature change for cooling
    t_cooling.comment.put("T cooling")
    t_cooling.expression.put("-(H-B)*D")

    # Temperature change for heating.
    # No heating if on_off is OFF.
    t_heating.comment.put("T heating")
    # Since transform values are floats, we evaluate the "boolean"
    # as True if value>0.5 (and False if <= 0.5).
    t_heating.expression.put("F>0.5? (A-B)*E: 0")

    # Temperature change for uniform random noise with amplitude G.
    # RNDM returns uniform random number between 0 .. 1.
    # Keep in mind that smoothing can contribute to overshoot since non-zero
    # smoothing introduces some lag in the computed temperature.
    t_noise.comment.put("T noise")
    t_noise.expression.put("G * (RNDM-0.5)*2")

    # smoothing fraction (0: only new values, 1: no new values)
    smoothing.comment.put("smoothing")
    smoothing.current_value.put(0.001)  # 0.0 .. 1.0

    # User changes the desired temperature here.
    # This value will be passed _to_ the epid record.
    setpoint.comment.put("setpoint")
    setpoint.current_value.put(starting)
    setpoint.output_pv.put(epid.final_value.pvname)  # _to_ epid

    # readback: the current temperature.
    # During steady-state with the heater on, T_heating should balance T_cooling
    # and the power_fraction should remain steady.  Any variation in
    # power_fraction should be in response to the effect of T_noise variations.
    readback.comment.put("readback")
    readback.expression.put("min(max(B, H+I+J+K), A)")

    # readback - setpoint
    difference.comment.put("difference")
    difference.expression.put("N-M")

    # Is the heater "at temperature"?
    at_temperature.comment.put("at temperature")
    at_temperature.expression.put("abs(O)<=C")

    # process epid after transform
    transform.forward_link.put(epid.prefix)

    # epid record wll be processed after transform by FLNK
    epid.scanning_rate.put("Passive")  # do not change

    # If Kp too low, controller will be slow to respond.
    # If Kp>=0.004, controller will oscillate!
    # If Ki is smaller, controller is slower to reach setpoint.
    # If Ki is larger, controller reacts to noise when at temperature
    # Kd is 0 for non-mechanical systems (do not change from zero)
    epid.proportional_gain.put(0.000_04)  # Kp
    epid.integral_gain.put(0.5)  # Ki
    epid.derivative_gain.put(0.0)  # Kd
    epid.high_limit.put(1.0)  # enforce 0.0 .. 1.0 range as power_fraction
    epid.low_limit.put(0.0)  # enforce 0.0 .. 1.0 range as power_fraction
    epid.low_operating_range.put(0)  # low == high: permissive
    epid.high_operating_range.put(0)  # low == high: permissive

    # This is readback: PID output is changed to minimize (readback-setpoint)
    epid.controlled_value_link.put(readback.current_value.pvname)

    # Since power_fraction.current_value _pulls_ epid.output_value,
    # do not configure epid.output_location to _push_ the value.
    epid.output_location.put("")  # leave this empty

    transform.scanning_rate.put(".1 second")

Make the changes:

[2]:

heater_simulator_setup("gp:userTran1", "gp:epid1")

heater as positioner#

The heater temperature may be described in ophyd as a positioner, enabling temperature scans and (thermal profiles)[https://en.wikipedia.org/wiki/Thermal_profiling] such as ramp, soak, cool.

Remember to turn on the power to the heater, or the controller will not come up to temperature.

[3]:

from apstools.devices import PVPositionerSoftDoneWithStop
from ophyd import Component
from ophyd import EpicsSignal

class HeaterPositioner(PVPositionerSoftDoneWithStop):
    """
    Simulated Heater as ophyd Positioner.

    PVs for setpoint (``.M``) and readback (``.N``) are defined through keyword
    arguments.  It is not necessary to create Components for them here. The
    ``PVPositionerSoftDoneWithStop`` will create components for both of these.

    EXAMPLE::

        temperature = HeaterPositioner(
            "gp:userTran1", name="temperature",
            setpoint_pv=".M", readback_pv=".N")
    """
    tolerance = Component(EpicsSignal, ".C", kind="config")
    on_off = Component(EpicsSignal, ".F", kind="config")

temperature = HeaterPositioner(
    "gp:userTran1", name="temperature",
    setpoint_pv=".M", readback_pv=".N")
temperature.wait_for_connection()
temperature.on_off.put(1)

Move the positioner using it’s ophyd controls:

[4]:

temperature.move(100)

[4]:

MoveStatus(done=True, pos=temperature, elapsed=20.2, success=True, settle_time=0.0)

[5]:

temperature.move(28.5)

[5]:

MoveStatus(done=True, pos=temperature, elapsed=21.3, success=True, settle_time=0.0)

Minimal setup of bluesky to demonstrate the positioner with the RunEngine. Use the plan_stubs to build a plan.

[6]:

from bluesky import RunEngine, plan_stubs as bps

RE = RunEngine()

Make a custom report() function.

[11]:

def report():
    sp = temperature.setpoint.get()
    rb = temperature.position
    print(
        f"  rb={rb:.3f}"
        f"  sp={sp:.3f}"
        f"  diff={rb-sp:.3f}"
    )

report()

  rb=28.604  sp=28.500  diff=0.104

Create a custom plan to demonstrate (with the RunEngine) the temperature object as a positioner. Use the move relative (bps.mvr()) plan stub to change the temperature.

[12]:

def my_plan(step=None):
    step = step or temperature.tolerance.get()

    report()

    yield from bps.mvr(temperature, step)
    report()

    yield from bps.mvr(temperature, -step)
    report()

First, run the plan, taking a 5 degree step. The step size is chosen to be well beyond the temperature.tolerance value.

[13]:

RE(my_plan(5))

  rb=28.595  sp=28.500  diff=0.095
  rb=32.750  sp=33.595  diff=-0.845
  rb=28.746  sp=27.750  diff=0.996

[13]:

()

There are problems with this plan:

setpoint did not change by the exact step value
setpoint did not return to the starting value as expected

That’s because bps.mvr() set the new setpoint as a relative change from the temperature’s position at the time a new step was requested.

Change the plan to control the temperature setpoint instead.

Before running the revsed plan, set the temperature back to 28.5.

[14]:

def my_plan(step=1):
    report()

    yield from bps.mvr(temperature.setpoint, step)
    report()

    yield from bps.mvr(temperature.setpoint, -step)
    report()

temperature.move(28.5)
RE(my_plan(5))

  rb=27.877  sp=28.500  diff=-0.623
  rb=27.877  sp=33.500  diff=-5.623
  rb=27.877  sp=28.500  diff=-0.623

[14]:

()

Good, the setpoint returned to the starting 28.5. New problems appeared:

The starting readback was not the expected 28.5
The readback did not change by the expected amount.

Both of these problems have the same cause, the temperature readback had not yet settled to the setpoint when the setpoint was next updated. While the readback was within the tolerance, the readback will be closer to the setpoint if we either:

wait a short time longer
reduce the tolerance

If the tolerance is too small, the temperature noise will cause the at_temperature value to fluctuate even if the setpoint is not changing.

Rather than reduce the tolerance, add the settling time into the plan (using bps.sleep()).

Also, switch back to absolute moves (bps.mv()) and controlling temperature (not temperature.setpoint).

[15]:

def my_plan(step=1, settle_time=5):
    report()
    starting = temperature.setpoint.get()

    yield from bps.mv(temperature, starting + step)
    yield from bps.sleep(settle_time)
    report()

    yield from bps.mv(temperature, starting)
    yield from bps.sleep(settle_time)
    report()

RE(my_plan(5))

  rb=28.326  sp=28.500  diff=-0.174
  rb=32.872  sp=33.500  diff=-0.628
  rb=28.934  sp=28.500  diff=0.434

[15]:

()

Good. Now we see how to change the temperature in a step-wise manner, using bps.mv() and a settling time using bps.sleep().

ramp the temperature#

When the tmperature is ramped, the setpoint is varied only some predetermined trajectory until some end condition is satisfied. One type of ramp is at constant rate where the temperature changes by a constant number of degrees / second (or minute).

In the ramp() plan below, the duration of the ramp is computed, given the start and final positions and the ramp rate. The setpoint is continually updated during the ramp (at the given period) to keep the temperature changing. The setpoint is a linear function of elapsed time, the start & final temperatures, and the direction of the ramp. The ramp ends when the duration has elapsed. One move to the final temperature finishes the ramp.

[22]:

import time

def ramp(positioner, final, rate=1, period=0.1, settle_time=5):
    start = positioner.setpoint.get()
    direction = 1 if final > start else -1
    duration = abs(final - start) / rate

    t0 = time.time()
    t_update = t0 + period
    while time.time() < t0 + duration:
        t = time.time()
        if t >= t_update:
            t_update += period
            value = start + direction * rate * (t - t0)
            yield from bps.mv(positioner.setpoint, value)
        yield from bps.sleep(period / 10)

    yield from bps.mv(positioner.setpoint, final)
    yield from bps.sleep(settle_time)

Ramp the temperature from the current value to 40.0.

[23]:

report()
RE(ramp(temperature, 40, settle_time=10))
report()

  rb=28.500  sp=28.500  diff=0.000
  rb=39.583  sp=40.000  diff=-0.417

There are new problems:

the readback temperature after the first ramp is not within tolerance (1.0)
during the ramp(), the following error was much greater than the temperature.tolerance and the controller was not at_temperature . The controller did not keep up with the setpoint changes

Either reduce the ramp rate or change the PID Kp and Ki terms to keep the following error smaller.

[25]:

# return to room temperature
temperature.move(28.5)
RE(bps.sleep(10))

# ramp again, slower this time
report()
RE(ramp(temperature, 40, rate=0.2, settle_time=10))
report()

  rb=28.752  sp=28.500  diff=0.252
  rb=39.607  sp=40.000  diff=-0.393

Good. The following error was smaller in this last ramp and the final readback was within tolerance of the setpoint.

Thermal profile#

With the ramp() plan above, a plan to ramp, soak, and cool is possible. Remember to control the ramp rate.

[26]:

def ramp_soak_cool(positioner, high, rate=1, soak_time=10):
    start = positioner.setpoint.get()
    report()

    print(f">>> ramp from {start:.2f} to {high:.2f}")
    yield from ramp(positioner, high, rate=rate)
    report()

    print(f">>> soak for {soak_time:.2f} s")
    yield from bps.sleep(soak_time)
    report()

    print(f">>> ramp from {high:.2f} to {start:.2f}")
    yield from ramp(positioner, start, rate=rate)
    report()

Run the new plan:

[27]:

RE(ramp_soak_cool(temperature, 55, rate=0.2, soak_time=15))

  rb=40.146  sp=40.000  diff=0.146
>>> ramp from 40.00 to 55.00
  rb=54.550  sp=55.000  diff=-0.450
>>> soak for 15.00 s
  rb=54.890  sp=55.000  diff=-0.110
>>> ramp from 55.00 to 40.00
  rb=40.526  sp=40.000  diff=0.526

[27]:

()