Implementing Euler and Tustin ODE Solvers in Siemens SCL

Overview

Ordinary differential equations (ODEs) describe the dynamic behavior of virtually every physical plant: motor speed, tank level, temperature loop, and current/voltage in an RLC network. In a Siemens PLC, the controller cannot integrate symbolically the way an engineer does with a Laplace transform on paper; the CPU must produce a discrete-time approximation of the continuous-time response. This article documents four production-tested SCL (Structured Control Language) function blocks that integrate first-order, second-order, and differentiator transfer functions in real time on S7-1200 and S7-1500 controllers running TIA Portal V16 or later.

Three blocks implement explicit Euler integration and one block implements the Tustin (bilinear) transform. Both techniques belong to the family of numerical methods for ordinary differential equations, the canonical reference for which is the Wiley text Numerical Methods for Ordinary Differential Equations. Each block was refined to set an explicit error flag whenever the user enters a coefficient that would cause division by zero, and to zero the integrators on the rising edge of a REST request.

Mathematical Background

All four function blocks solve variants of a single canonical form. The first-order plant

Y(s)/F(s) = 1 / (A·s + B)

is the generic form of a low-pass filter. The casual differentiator

Y(s)/F(s) = C·s / (A·s + B)

is the D-term of an ideal PID controller written in parallel form. The second-order system

Y(s)/F(s) = 1 / (A·s² + B·s + C)

describes a lightly damped mechanical resonance, an LC filter, or a positioning axis. By scaling A, B, and C, the same code covers critically damped, under-damped, and over-damped responses.

Euler Integration

Explicit (forward) Euler replaces the differential operator with the finite difference dy/dt ≈ (Y(k+1) − Y(k)) / T, giving the recurrence

Y(k+1) = Y(k) + T · F(k, Y(k))

with local truncation error O(T²) and global error O(T). It is conditionally stable: for the first-order plant with time constant τ, the block remains stable only if the sample time T < 2·τ. Engineers using Euler in production must therefore halve the controller's task time whenever they double the dominant time constant.

Tustin (Bilinear) Transform

The Tustin transform substitutes

s ≈ (2/T) · (z − 1) / (z + 1)

which maps the left-half s-plane onto the interior of the unit circle and yields unconditional stability for any positive T. For the first-order plant this gives the difference equation

(2B + C·T)·Y(k+1) = A·T·F(k+1) + A·T·F(k) + (C·T − 2B)·Y(k)

The mapping is exact at DC and is exact at half the Nyquist frequency, with smooth amplitude roll-off in between. Engineers typically use Tustin for plant models that will be simulated in the same scan as the controller, where stability over a wide T range outweighs the small extra arithmetic cost.

SCL Environment and Prerequisites

The blocks were developed against the SCL grammar that ships with:

• STEP 7 V5.6 SP2 (S7-300/400 target)
• STEP 7 in TIA Portal V16, V17, V18, V19 (S7-1200/1500 target)

The SCL language manual for S7-1500 is available on the Siemens support portal as "S7-1500/ET 200MP - Programming and Operating Manual - SCL"; the equivalent for S7-1200 is "S7-1200 Programmable Controller - System Manual". Both describe the TIME data type, the TIME_TO_DINT and DINT_TO_REAL conversion functions, and the REAL arithmetic library used by every block below.

First-Order System via Euler (FB1)

FB1 solves Y/F = 1 / (A·s + B). The integrator is updated every INTERVAL, defaults to 1 s, and may be re-tasked down to 10 ms on S7-1500 and 100 ms on S7-1200. The block uses YOLD as the static integrator state, applies an explicit Euler step, then writes the new value back into the static for the next scan.

FUNCTION_BLOCK FB1
TITLE = 'ODE_EULER'
// Y      1
// -- = ----------
// F    A*S + B
VERSION: '2.0'
AUTHOR:  HD
NAME:    first
FAMILY:  FORUM_E
VAR_INPUT
  F        : REAL;
  A        : REAL;
  B        : REAL;
  INTERVAL : TIME  := T#1S;
END_VAR
VAR_IN_OUT
  REST     : BOOL  := FALSE;
END_VAR
VAR_OUTPUT
  Y        : REAL;
  error    : BOOL;
END_VAR
VAR
  YNEW     : REAL  := 0.0;
  YOLD     : REAL  := 0.0;
  K        : REAL  := 0.0;
  DELTA    : REAL  := 0.0;
  T_INTERNAL: REAL := 1.0;
END_VAR
BEGIN
  T_INTERNAL := DINT_TO_REAL(TIME_TO_DINT(INTERVAL)) / 1000.0;
  IF REST = 1 THEN
    YOLD  := 0;
    YNEW  := 0;
    Y     := 0;
    OK    := FALSE;
    error := 1;
  ELSIF A = 0 THEN
    YOLD  := 0;
    YNEW  := 0;
    Y     := 0;
    OK    := FALSE;
    error := 1;
  ELSE
    DELTA := (F - (B * YOLD)) / A;
    YNEW  := YOLD + (DELTA * T_INTERNAL);
    YOLD  := YNEW;
    Y     := YNEW;
    OK    := TRUE;
    error := 0;
  END_IF;
END_FUNCTION_BLOCK

Calling Convention

1. Instantiate the FB in a cyclic OB with a unique instance DB.
2. Drive F with the plant input (REAL in engineering units).
3. Set A to the time-constant τ and B to the static gain denominator coefficient.
4. Set INTERVAL equal to the OB cycle or the desired update period; the conversion TIME_TO_DINT(...)/1000.0 produces seconds in REAL.
5. Pulse REST for one cycle to clear the integrator; the block holds error=1 for that cycle.

Casual Differentiator via Euler (FB2)

FB2 implements Y/F = C·s / (A·s + B), the derivative branch of an ideal PID controller. The same time constant A/B filters the raw derivative to keep it realizable, and C sets the high-frequency gain. Because the block contains an explicit derivative on F, it must use FOLD to approximate dF/dt over the same interval.

FUNCTION_BLOCK FB2
TITLE = 'ODE_EULER'
// Y   C*S
// -- = ------
// F   A*S + B
VERSION: '2.0'
AUTHOR: HD
NAME:   DIFF
FAMILY: FORUM_E
VAR_INPUT
  F        : REAL;
  A        : REAL  := 1.0;
  B        : REAL  := 1.0;
  C        : REAL  := 1.0;
  INTERVAL : TIME  := T#1S;
END_VAR
VAR_IN_OUT
  REST     : BOOL  := FALSE;
END_VAR
VAR_OUTPUT
  Y        : REAL;
  error    : BOOL;
END_VAR
VAR
  FOLD     : REAL  := 0.0;
  YNEW     : REAL  := 0.0;
  YOLD     : REAL  := 0.0;
  DELTA    : REAL  := 0.0;
  T_INTERNAL: REAL := 1.0;
END_VAR
BEGIN
  T_INTERNAL := DINT_TO_REAL(TIME_TO_DINT(INTERVAL)) / 1000.0;
  IF REST = 1 THEN
    YOLD  := 0;
    YNEW  := 0;
    Y     := 0;
    OK    := FALSE;
    error := 1;
  ELSIF A = 0 THEN
    YOLD  := 0;
    YNEW  := 0;
    Y     := 0;
    OK    := FALSE;
    error := 1;
  ELSE
    DELTA := ((C * (F - FOLD)) / T_INTERNAL - (B * YOLD)) / A;
    YNEW  := YOLD + (DELTA * T_INTERNAL);
    YOLD  := YNEW;
    FOLD  := F;
    Y     := YNEW;
    OK    := TRUE;
    error := 0;
  END_IF;
END_FUNCTION_BLOCK

Second-Order System via Euler (FB3)

FB3 implements Y/F = 1 / (A·s² + B·s + C) by writing the system as two coupled first-order states X1 = Y and X2 = Ẏ, then integrating each with explicit Euler. This is the standard phase-variable decomposition used in control theory textbooks.

FUNCTION_BLOCK FB3
TITLE = 'ODE_EULER'
// Y       1
// -- = -------------
// F    A*S^2 + B*S + C
VERSION: '2.0'
AUTHOR: HD
NAME:   SECOND
FAMILY: FORUM_E
VAR_INPUT
  F        : REAL;
  A        : REAL  := 1.0;
  B        : REAL  := 1.0;
  C        : REAL  := 1.0;
  INTERVAL : TIME  := T#1S;
END_VAR
VAR_IN_OUT
  REST     : BOOL  := FALSE;
END_VAR
VAR_OUTPUT
  Y        : REAL;
  error    : BOOL;
END_VAR
VAR
  X1       : REAL  := 0.0;
  X2       : REAL  := 0.0;
  X1OLD    : REAL  := 0.0;
  X2OLD    : REAL  := 0.0;
  DELTA1   : REAL  := 0.0;
  DELTA2   : REAL  := 0.0;
  T_INTERNAL: REAL := 1.0;
END_VAR
BEGIN
  T_INTERNAL := DINT_TO_REAL(TIME_TO_DINT(INTERVAL)) / 1000.0;
  IF REST = 1 THEN
    X1     := 0.0;
    X2     := 0.0;
    X1OLD  := 0.0;
    X2OLD  := 0.0;
    Y      := 0;
    OK     := FALSE;
    error  := 1;
  ELSIF A = 0 THEN
    X1     := 0.0;
    X2     := 0.0;
    X1OLD  := 0.0;
    X2OLD  := 0.0;
    Y      := 0;
    OK     := FALSE;
    error  := 1;
  ELSE
    DELTA2 := (F - (C * X1OLD) - (B * X2OLD)) / A;
    X2     := X2OLD + (DELTA2 * T_INTERNAL);
    DELTA1 := X2;
    X1     := X1OLD + (DELTA1 * T_INTERNAL);
    X2OLD  := X2;
    X1OLD  := X1;
    Y      := X1;
    OK     := TRUE;
    error  := 0;
  END_IF;
END_FUNCTION_BLOCK

Stability Caveat

For a second-order plant A·s² + B·s + C, the explicit Euler update matrix has eigenvalues that leave the unit circle when T > 2·A/B. Always satisfy T < A·ωn·2 / B (where ωn = sqrt(C/A)) and run the block from a fast OB (OB35 at 100 ms or OB1 at 10 ms on S7-1500) when the natural frequency approaches 1 Hz.

First-Order Process via Tustin (FOP)

The Tustin implementation avoids the conditional-stability problem of Euler at the cost of a small amount of algebra at every cycle. The block stores the previous input U0 and output Y0, then computes Y1 from the bilinear form of the transfer function Y/U = A / (B·s + C).

FUNCTION_BLOCK FOP
TITLE = 'First Order Process'
// Y         A
// -- = ----------
// U       B*s + C
//
// Use Tustin method to approximate
//       2     z-1
// s = ---*----------
//       T     z+1
//
//  Y               A*T*z + A*T
//  -- = ----------------------------------
//  U          (2B + C*T)*z + (C*T - 2B)
//
// (2B + C*T)*y(k+1) = A*T*u(k+1) + A*T*u(k) + (C*T - 2B)*y(k)
//
VERSION: '0.0'
AUTHOR: ZLiang
NAME:   FOP
FAMILY: PROCESS
VAR_INPUT
  I_U   : REAL;
  I_A   : REAL;
  I_B   : REAL;
  I_C   : REAL  := 1.0;
  I_Ts  : TIME  := T#100MS;
END_VAR
VAR_IN_OUT
  IQ_RST: BOOL  := FALSE;
END_VAR
VAR_OUTPUT
  Q_Y   : REAL;
  Q_ERR : BOOL;
END_VAR
VAR
  Y0  : REAL := 0.0;
  Y1  : REAL := 0.0;
  U0  : REAL := 0.0;
  Ts  : REAL;
END_VAR
BEGIN
  Ts := DINT_TO_REAL(TIME_TO_DINT(I_Ts)) / 1000.0;
  IF IQ_RST THEN
    Y0    := 0.0;
    Y1    := 0.0;
    U0    := 0.0;
    Q_ERR := FALSE;
    IQ_RST:= FALSE;
  ELSIF (I_B < 0.005) AND (I_B > -0.005) THEN
    Q_ERR := TRUE;
    Y0    := 0.0;
    Y1    := 0.0;
  ELSE
    Y1    := ((I_A * Ts * I_U) + (I_A * Ts * U0) +
              ((2.0 * I_B) - (I_C * Ts)) * Y0) /
             ((2.0 * I_B) + (I_C * Ts));
    U0    := I_U;
    Y0    := Y1;
    Q_ERR := FALSE;
  END_IF;
  Q_Y := Y1;
END_FUNCTION_BLOCK

Why Tustin over Euler for Plant Models

The Tustin form matches the plant exactly at DC, so the steady-state gain of the simulated plant is always 1.0 (for A/C = 1) regardless of sample time. When the same scan computes the controller and the plant, Tustin avoids the drift that explicit Euler introduces into long-running simulations, particularly for integrator-like plants with B ≈ 0 (the block guards that case with the |I_B| < 0.005 check and sets Q_ERR).

Method Comparison

Stability and Accuracy Analysis

For the first-order plant A·s + B, explicit Euler produces the closed-form iteration

Y(k+1) = (1 − B·T/A)·Y(k) + (T/A)·F(k)

which is bounded-input bounded-output (BIBO) stable iff |1 − B·T/A| < 1, i.e. T < 2·A/B. Engineers tuning a PID loop that includes this block should keep the cycle time below 50 % of the dominant plant time constant to leave a comfortable stability margin.

The Tustin block's iteration is

Y(k+1) = [(C·T − 2B)/(2B + C·T)]·Y(k) + [A·T/(2B + C·T)]·(F(k+1) + F(k))

with pole at (2B − C·T)/(2B + C·T), which is strictly inside the unit circle for all positive B, C, T. Amplitude error at ω = ωn is 0 dB; phase error at ω = ωn is below 5° for T < 0.2·τ.

Commissioning and Verification

1. Compile and download the four FBs into the PLC project. Each FB creates its own instance DB on the first call.
2. Open the instance DB online and confirm the YOLD, X1OLD, X2OLD, FOLD, Y0, and U0 statics initialize to 0.0. Any non-zero starting condition produces a step transient on the first scan.
3. Apply a step on F and trace Y with a 1 kHz sampling trace (S7-1500 supports this in the Online & Diagnostics view). Verify the steady-state value matches the analytical Y(∞) = F·A/B (first order) or Y(∞) = F·A/C (second order) within 1 %.
4. Halve the OB period and confirm the simulated step response is unchanged. With Euler, halving the period reduces the error by 2×; with Tustin, the response is essentially identical to four significant figures.
5. Force A = 0 on FB1, FB2, or FB3 and confirm the block asserts error = 1 and clamps Y = 0 without throwing an OB priority-class exception. The PLC remains in RUN even with a forced zero coefficient.
6. Force I_B near 0 on FOP and confirm Q_ERR is raised. Reset IQ_RST to clear.

Integration Patterns in a TIA Portal Project

Engineers typically wire the four blocks into a CFC chart or a sequence of SCL calls in OB35. A recommended layout for a position-loop simulation:

1. FB1 simulates the first-order current loop, with A = τ_i, B = 1, and F driven by the controller output.
2. FB3 simulates the second-order mechanical plant, with A = J, B = B_friction, C = K_spring, and F = current-loop output.
3. FB2 forms the D-term of a parallel PID, with A = τ_d + 0.05 to clamp derivative noise and C = K_d.
4. FOP replicates the same first-order plant using Tustin, sitting in parallel with FB1 to validate the simulation against the analytical response during commissioning.

Error Handling Matrix

Frequently Asked Questions