# Simulating drift, stiction and lost motion

*Read part 1 of this article here.*

Adding realism in the form of random events to an otherwise deterministic, dynamic process simulator makes simulation more realistic for evaluating control structures. The following is an excerpt from the ISA-published book, *Nonlinear Model-based Control: Using First-principles Models in Process Control*.

**Get your subscription to Control's tri-weekly newsletter.**

## Simulating ever-changing drifts

The simulation should include mechanisms that generate ever-changing, environmental influences, input conditions, calibration drifts and process attributes. There are four categories of such variables:

- Environmental influences include ambient heat losses due to weather or relative humidity in air intakes;
- Input conditions include the BTU content of a fuel, raw material composition or flow rates;
- Calibration drifts can be on any process measurement, and related to instrument warmup, instrument aging and sensor fouling, etc; and
- Process attributes, including distillation tray efficiency, orifice erosion, screen plugging, heat transfer area fouling and catalyst reactivity.

Though many control theories identify “a disturbance” as one additive bias to a process output, there are many disturbances that don’t act as a measurement bias to the simulator output. They should be modeled as process-input influences.

Certainly, it’s easy to have a computer make step or ramp changes in an input, and some influences do switch on and off. However, this is not commonly the way nature devises it. More consistent with reality, the input pattern should be characterized by a random walk about a mean value—a stochastic pattern. Let’s look at how to include disturbances as auto-regressive, moving- average (ARMA) simulator inputs.

To create a simulator, first define the equation or procedure to obtain y_{true}, the truth about nature, from x_{desired}, the desired value of the influence:

y_{true} = f (x_{desired},p_{ideal},e_{ideal},c_{ideal},in_{ideal}) (1)

Though shown as an equation, the process simulator will likely be a complicated set of functions and subroutines. The x_{desired} represents controlled and other measured variables. The p_{ideal} represents ideal values of the process or equipment attributes (catalyst reactivity, tire air pressure, heat exchanger fouling, fluid flow pressure loss coefficient, etc). The e_{ideal} represents ideal values of the environmental influences (temperature, %RH, etc). The c_{ideal} represents the ideal values of the instrument calibrations (for instance, 4-20 mA generates 3 - 15 psia, which makes a valve stem range from 0 - 100% of stroke). The in_{ideal} represents the ideal values of the process inflows (temperature, raw material composition, etc).

Before using the model, first perturb each variable in each of the three categories with appropriate disturbances (drifts, d), and, if sensible, add noise, n. Then use the model to calculate the true response. Finally, perturb the response with measurement drift and noise:

x_{actual} = x_{desired} + d_{x} + n_{x} (2a)

p_{actual} = p_{ideal} + d_{p} + n_{p} (2b)

e_{actual} = e_{ideal} + d_{e} + n_{e} (2c)

c_{actual} = c_{ideal} + d_{c} + n_{c} (2d)

in_{actual} = in_{ideal} + d_{in} + n_{in} (2e)

y_{model} = f (x_{actual},p_{actual},e_{actual},c_{actual},in_{actual}) (2f)

y_{measured} = y_{model} + d_{y} + n_{y} (2g)

Although the disturbances, d, and noise, n, are added, it might be more appropriate to make them a multiplicative factor, when you anticipate that the noise level or perturbation magnitude scales with the input value. Although these drift and noise values change with each sampling, for clarity, the subscript counter, i, is not included.

Looking at measurement sequence or time, there are two basic concepts of variation. One is independent variation, in which sample-to-sample variation is uncorrelated. This is often termed noise.

With noise, what causes one perturbation from ideal doesn't persist, and the next perturbation is the result of independent causes. The other mechanism for variation reflects a cause that has persistence, in which sequential perturbations are correlated. The ever-changing disturbances, labeled as d in the previous equations, represent auto-correlated, time-dependent trends, or drifts, in an influence or model coefficient. The symbol n will represent uncorrelated and independent noise.

Many influences on experiments have persistence. For example, if a cloud covers the sun, the shadow persists for a while. If the blocked sunlight influences a temperature disturbance, then a negative d-value at one sampling (in the shadow of the cloud) generates a negative d-value at the subsequent sampling (as the shadow persists).

For dynamic, time-dependent simulators, a common method to simulate how disturbances (bias or systematic error in measurements, as well as uncontrolled influences on the process) change over time is to considering that they’re driven by random events that have persistence.

A drifting influence on the j-th variable with a first-order persistence driven by NID(0,σ) noise is modeled as:

dd_{j}(t) + d_{j}(t) = n(t) = σ_{d} √(-2ln (r_{(1,i)})) sin (2πr_{(2,i)}) (3)

dt

Again, the subscript i represents the sample number—the time interval. Expanding this differential with a forward-finite difference and rearranging leads to the ARMA model for d_{j}(t):

d_{(j,i)} = (1 – λ) d_{(j,i-1)} + λσ_{d} √(-2ln(r_{(1,i)})) sin (2πr_{(2,i)}) (4)

Where,

λ = 1 – e^{(-∆t/τ)} (5)

and ∆t is the simulation time step.

If ∆t < τ, which is usually the case in simulation, then the simplified approximation λ = ∆t/τ can be used.

When creating a simulation, the user chooses a time-constant for the persistence that’s reasonable for the effects considered, and a σ-value that would make the disturbance have a reasonable variability. Since the first-order persistence tempers the response to the random influences, the driver, σ_{d}, must be greater than the desired variation on the disturbance.

When creating a simulation with a stochastic influence, the user chooses a time-constant for the persistence that’s reasonable for the effects considered, and a σ_{d}-value that would make the disturbance have a reasonable variability. At each sampling, Equation 4 provides the perturbation value for the variables x, p, e, c or in that are used in all of Equation 2. In its first use, the values for d_{j} in Equation 4 should be initialized with zero.

How should one choose values for λ and σ_{d}? First, consider the time-constant, τ, in Equation 3. It represents the time-constant of the persistence of a particular influence. It's roughly τ ≈ 1/3 of the lifetime of a persisting event because the solution to the first-order differential equation indicates that, after three time-constants, d_{j} finished 95% of its change toward the final value. If you consider that the shadow of a cloud persists for six minutes, then the time-constant value is about two minutes. Or you could use τ ≈ 1/4 of the persistence period representing 98% of the time to change. Once you choose a τ-value that matches your experience with nature and decided on a time interval for the numerical simulation, ∆t, calculate λ from Equation 5.

To determine the value for σ_{d}, propagate variance in Equation 4. You don’t have to do it. The result is Equation 6. Use your choice of σ_{z} (and λ, which is dependent on your choices for ∆t and τ) to determine the value for σ_{d}. The subscript z represents any of the variables x, p, e, c or in:.

σ_{d} = σ_{z }√((2 – λ) /λ) (6)

To choose a value for σ_{z}, the resulting variability on the z-variable, consider a fluctuation range for the disturbance. You must have a feel for what’s reasonable to expect for the situation that you’re simulating. For instance, if it’s barometric pressure, the normal local range of low-to-high might be 29 - 31 inches of mercury (inHg). If it’s outside temperature in the summer, it might be 70 - 95 oF. If it’s catalyst activity coefficient, it might be from 0.50 to 0.85. The disturbance value is expected to wander within those extremes. Using the range, R, as:

R = HIGH – LOW (7)

And the standard deviation, σ, as approximately one-fifth of the range, then:

σ_{d} = (R/5) √((2-λ) / λ) (8)

To generate stochastic inputs for dynamic simulators, choose a time-constant for persistence of events and a range of the disturbance variable. Use Equation 5 to calculate λ, then Equations 6 and 8 to calculate σ_{d}. At each simulation time interval, use Equation 4 to determine the perturbation to a variable.

Figure 1 is an illustration of one 400-minute realization of a stochastic variable calculated as above. The time-constant is 40 minutes, the nominal input value is 20, and the range is three units. In this illustration, notice that the variable averages about its nominal value of 20, and that the difference between the high and low values is nearly three. During a period between 150 and 275 minutes, the value is below the average, indicating a persistence of 275 - 150 = 125 minutes, which is about three times the 40-minute time-constant. Some persistence values are shorter, while some are longer.

This approach modeled a drifting influence as a first-order response to a Gaussian distributed driver. You can certainly can model drifts as higher-order or driven by alternate, forcing-function distributions. However, my view is that this adequately mimics how I consider randomly driven drifts in natural disturbances to respond, and it provides a relatively simple method to execute and parameterize.

## Simulating stiction

Stiction is a slip-stick phenomenon common to pneumatically activated flow-control valves. The valve stem must slide through the valve body to open and close the inner valve. This hole in the body permits process material to leak out, so there’s packing material between the stem and body. If the packing is too loose, process material can slip out. But if it’s too tight, the packing can grab and hold the valve stem.

Normally, when the controller wants to open the a valve, the air pressure on the actuator changes, and the force imbalance between the actuator and the spring permits the stem to move to a position that makes the forces balance. But, if the packing grabs the stem, the stem will not move until the force imbalance overcomes the static friction. Then, the valve stem starts moving. And, since sliding friction is less than static friction, the stem jumps to a new position, stopping at a random position between the current location and the ideal target, where the static friction again overcomes the actuator-to-spring force imbalance.

A convenient way to simulate stiction is to separately model the target stem position (that desired by the controller), x_{target}, and the actual stem position, x_{actual}. If the difference between the target and actual is less than a threshold difference, then the actual position doesn't change. However, if the difference exceeds the threshold, ∆, then the actual jumps to a random spot near to the target:

IF (| x_{target} – x_{actual} | < ∆ ) THEN

(x_{actual}: = x_{actual} ) (9)

IF (| x_{t} – x_{a} | ≥ ∆) THEN

(x_{a}:= x_{a} + (1+r) (x_{t} – x_{a}) / 2) (10)

In Equation 10, the target and actual subscripts are replaced with “t” and “a.” The r represents a uniformly distributed random number in the 0 to 1 interval. The (1 + r) / 2 functionality in this model permits the valve stem to jump at least halfway to the target prior to stopping, randomly. The := notation indicates a computer code assignment statement, in which the “prior” and “new” subscripts are unnecessary.

The user specifies the threshold, ∆. Certainly, the user can make the ∆-value dependent on the stem position (modeling the local impact of defects) or the jump-to range larger or smaller.

Ideally, in a single-input, single-output (SISO) application, stiction is shown by a square wave (alternating steps) in the CV and a corresponding sawtooth wave (alternating ramps) in the MV. The pattern might not be exactly regular, and may be shifted by a delay or tempered by a filter. The evidence may be more complicated in a multivariable process and with continually drifting disturbance effects.

## Simulating lost motion

Lost motion happens when there’s play in mechanical linkages, which can be on a rotary or butterfly valve in the process industry, or in positioning devices in mechanical or robotic devices. With play or slop in linkage, the final element follows the target position—always behind by a small value. When the target position reverses, the final element doesn't change until the linkage slop is removed, and the final element can begin to move. Often, this is termed “deadband,” but lost motion is more appropriate because the term deadband normally refers to the band around a setpoint used to trigger on/off control action.

A convenient way to simulate this is to separately model the target, final-element position (that desired by the controller), x_{target}, and the actual position, x_{actual}. If the difference between the target and actual is less than a threshold difference, then the actual position doesn’t change. However, if the difference exceeds the threshold, ∆, then the actual

follows the target with a difference of ∆:

IF (|x_{target} – x_{actual} | < ∆)

THEN (x_{actual}:= x_{actual} ) (11)

IF (x_{t} – x_{a} ≥ ∆ )

THEN (x_{a}:= x_{t} – ∆) (12)

IF (x_{t }– x_{a} ≤ ∆)

THEN( x_{a}:= x_{t} + ∆) (13)

Don’t learn by memorizing these messages. Learn by implementing the techniques in a simulator and exploring their impacts.