Towards an Operational Semantics of the Simulation Engine of - - PowerPoint PPT Presentation

Towards an Operational Semantics of the Simulation Engine of Simulink

Alexandre Chapoutot

joint work with Olivier Bouissou (CEA LIST, LMeASI)

ENSTA Paristech

SYNCHRON 2011 December 1st

Motivation

2 / 22

Simulink and its industrial uses

Context

Simulink is a de facto standard in the industry to design embedded safety critical applications. e.g. it is used to design almost all the drive-by-wire systems.

Components of a drive-by-wire systems

A physical mechanism we want to control, e.g. a suspension. A software which implement the controller. ⇒ we have to deal with hybrid systems.

Current validation method

The main method to validate the design of embedded systems is based

• n simulation activity.

3 / 22

Classical design and validation methodology

m¨ x + c ˙ x + kx = u Mathematical model

zb relative_position 1 Integrator1 1 s Integrator 1 s Gain1 1/m Gain k c(t) 2 road 1 dot_zb ddot_zb Road_profile road Quarter_car road damping_force relative_position Controller relative_position damping_force

Controller definition

Identification Computer description experimental studies Theoretical and

Remarks

Mathematical models are an approximated descriptions of physical systems. The computer descriptions are studied with numerical tools.

4 / 22

Classical design and validation methodology

MIL: Model in the loop

Definition of a mathematical model of the plant and the controller.

Simulation platform

Test process Controller: model Plant Model Sensor model Actuator model Test vectors Expected output

Aim of the testing activity

Does the controller fulfil the specification? The set of tests will be used as “an oracle” in the next steps.

4 / 22

Classical design and validation methodology

SIL: Software in the loop

Implementation of the controller in a target language.

Simulation platform

Test process Controller: source code Plant Model Sensor model Actuator model Test vectors Expected output

Aim of the testing activity

Do the hand-written or generated code still fulfil the specification?

4 / 22

Classical design and validation methodology

PIL: Processor in the loop

Compilation of the controller and execute it on a virtual processor.

Simulation platform

Test process Controller: object code Virtual targeted CPU Plant Model Sensor model Actuator model Test vectors Expected output

Aim of the testing activity

Does the low-level implementation still fulfil the specification?

4 / 22

Considering Simulink as a language

Simulink: a graphical language

It allows to describe, as block-diagrams, a mixing of: Ordinary differential equations (ODE). Finite difference equations: single-rate or multi-rate period. Conditional executed equations: enabled or triggered.

Simulink: several parametrized semantics

the definition of a numerical solver used to solve ODEs; the detection of events: zero-crossing. So, it is a numerical approximation of the mathematical behavior.

What we should keep in mind!

The designer only refers to the Simulink’s semantics (i.e. graphical

• utput) to validate the development steps!

5 / 22

Formal verification of Simulink

Fact: simulations will never be complete and cannot bring strong guarantee on the designed software. Main question: How can we help designer to be more confident in its designed systems with formal methods?

Formal verification of Simulink: Our Goal

To compute invariant properties on the software taking into account a model of the plant. Our approach: abstract interpretation-based static analysis.

Formal verification of Simulink: Constraint

Formal verification methods should be adaptable for the design process in order to be more easily adopted by the end-user.

Formal verification of Simulink: Challenge

Understand the simulation engine of Simulink.

6 / 22

Simulink Language

7 / 22

A simple example: mathematical models

A semi-active suspension of a quarter-car model

Mass (m) Sensor k Controller zb − zr c(t) zb zr m = 250 kg k = 20000 N/m cmax = 16000 N/m/s cmin = 0

Mathematical model of the mechanical system

¨ zb = − 1 m

• k(zb − zr) + c(t)
• .

Mathematical model of the controller

c(t) =

• −cmax(zb − zr)

if (zb − zr)(˙ zb − ˙ zr) < 0 cmin if (zb − zr)(˙ zb − ˙ zr) ≥ 0 .

8 / 22

A simple example: Simulink implementation – 1

Mathematical model of the mechanical system

¨ zb = − 1 m

• k(zb − zr) + c(t)
• .

zb relative_position 1 Integrator1 1 s Integrator 1 s Gain1 1/m Gain k c(t) 2 road 1 dot_zb ddot_zb

Integrator block

Associated to a first order dynamic system:

• ˙

x(t) = input(t)

• utput(t) = x(t)

with x(0) = x0

9 / 22

A simple example: Simulink implementation – 2

Mathematical model of the controller

c(t) =

• −cmax(zb − zr)

if (zb − zr)(˙ zb − ˙ zr) < 0 cmin if (zb − zr)(˙ zb − ˙ zr) ≥ 0 . Implementation: (˙ zb − ˙ zr) is given by differentiating the sensor output.

damping_force 1 c_min c_max c Unit Delay z 1 Switch >= 0 Subtract Product1 Product Gain K relative_position 1

Discrete differentiation at rate 1/40 sec.

Closed-loop system

Connect the output of the plant to the input of the controller. Connect the output of the controller to the input of the plant.

9 / 22

A simple example: simulation results

Simulation settings

Duration of the simulation: 5 seconds. Variable step-size solver: ode23. Zero-crossing detection activated (non adaptive version).

2 4 5 · 10−2 0.1 Road profile 2 4 −0.1 − 5 · 10−2 5 · 10−2 Relative position

Remarks

The on-off controller make the suspension stable.

10 / 22

A simple example: simulation results

2 4 −1.5 −1 −0.5 ·104 Damping coefficient c(t) 100 200 2 4 Time vs. simulation loop iteration

Remark: numerical precision

The damping force oscillated because of tiny variations of the suspension. Consequence: finite precision may induce unexpected behaviors.

Remark: time

The time evolution is not homogeneous: Consequence: blocks depending on time may produce more or less

• utput following the time evolution e.g. sine block.

10 / 22

Operational Semantics

11 / 22

Simulink as an equation-based programming language

Input/Output relation

Each block defines the time invariant relation between its input and its

• utput.

Library Blocks Representation Equations Sources Input

1 In1

ℓ1

ℓ1 = in1 Constant

c Constant

ℓ1

ℓ1 = c Arithmetic Add

Add

ℓ1 ℓ2 ℓ3

ℓ3 = ℓ1 + ℓ2 Signal routing Switch

Switch ℓ1 ℓ2 ℓ3 ℓ4

ℓ4 = if (pr(ℓ2), ℓ1, ℓ3) Continuous-time Integrator

1/s Integrator

ℓ1 ℓ2

{ℓ2 = x; ˙ x = ℓ1; x(0) = init} Discrete-time Unit Delay

z 1 Unit Delay

ℓ2 ℓ1

{ℓ2 = d; ¯ d =S ℓ1; d(0) = init}

¯ d =S ℓ stands for at each t = kS, d = ℓ else it keeps its previous value.

12 / 22

Simulink as an equation-based programming language

Input/Output relation

Each block defines the time invariant relation between its input and its

• utput.

A core language of equations

e ::= r | ℓ | x | d | e1 ⋄ e2 | e1 ⊲ ⊳ e2 | if e1, e2, e3

• (1)

eq ::= ℓ :=S e | ℓ := e | ˙ x := e | ¯ d :=S e (2) p ::= eq | eq; p (3)

with r ∈ R, constant values; ℓ ∈ V, variables associated to a block output; x ∈ V, variables associated to continuous-time states; d ∈ V, variables associated to discrete-time states; ⋄ ∈ {+, −, ×, ÷}, arithmetic operations; ⊲ ⊳∈ {<, ≤, >, ≥, =, <>}, relational operations; S is the set of all the sampling times.

12 / 22

Simulink as an equation-based programming language

zb relative_position 1 Integrator1 1 s Integrator 1 s Gain1 1/m Gain k c(t) 2 road 1 dot_zb ddot_zb

1 2 3 12 13 14 15 x0 x1

damping_force 1 c_min c_max c Unit Delay z 1 Switch >= 0 Subtract Product1 Product Gain K relative_position 1

In red: evaluation

• rder of blocks

4 5 6 7 8 9 10 11 d For each block in the evaluation order, a simple translation gives:

ℓ1 = input ˙ x1 = ℓ15 ℓ2 = x1 ℓ3 = ℓ2 − ℓ1 ¯ d =S ℓ3 ℓ4 = d ℓ5 = ℓ3 − ℓ4 ℓ6 = 1/40 × ℓ5 ℓ7 = ℓ3 × ℓ6 ℓ8 = 0 ℓ9 = −16000 ℓ10 = if ℓ7 ≥ 0 then ℓ8 else ℓ9 ℓ11 = ℓ3 × ℓ10 ℓ12 = 20000 × ℓ3 ℓ13 = −ℓ12 − ℓ11 ℓ14 = 1/250 × ℓ13 ˙ x0 = ℓ14 ℓ15 = x0

12 / 22

Simulink as an equation-based programming language

Kinds of equations

A Simulink model is made of four kinds of functions: the output function; the update function of discrete-time states; the update function of continuous-time states;

ℓ1 = input ˙ x1 = ℓ15 ℓ2 = x1 ℓ3 = ℓ2 − ℓ1 ¯ d =S ℓ3 ℓ4 = d ℓ5 = ℓ3 − ℓ4 ℓ6 = 1/40 × ℓ5 ℓ7 = ℓ3 × ℓ6 ℓ8 = 0 ℓ9 = −16000 ℓ10 = if ℓ7 ≥ 0 then ℓ8 else ℓ9 ℓ11 = ℓ3 × ℓ10 ℓ12 = 20000 × ℓ3 ℓ13 = −ℓ12 − ℓ11 ℓ14 = 1 250 × ℓ13 ˙ x0 = ℓ14 ℓ15 = x0

Remark

These functions allow a state-space representation of a model.

12 / 22

Simulink as an equation-based programming language

Kinds of equations

A Simulink model is made of four kinds of functions: the output function; the update function of discrete-time states; the update function of continuous-time states;

ℓ1 = input ˙ x1 = ℓ15 ℓ2 = x1 ℓ3 = ℓ2 − ℓ1 ¯ d =S ℓ3 ℓ4 = d ℓ5 = ℓ3 − ℓ4 ℓ6 = 1/40 × ℓ5 ℓ7 = ℓ3 × ℓ6 ℓ8 = 0 ℓ9 = −16000 ℓ10 = if ℓ7 ≥ 0 then ℓ8 else ℓ9 ℓ11 = ℓ3 × ℓ10 ℓ12 = 20000 × ℓ3 ℓ13 = −ℓ12 − ℓ11 ℓ14 = 1 250 × ℓ13 ˙ x0 = ℓ14 ℓ15 = x0

Notation

We denote this function by g(t, x, d) .

12 / 22

Simulink as an equation-based programming language

Kinds of equations

A Simulink model is made of four kinds of functions: the output function; the update function of discrete-time states; the update function of continuous-time states;

ℓ1 = input ˙ x1 = ℓ15 ℓ2 = x1 ℓ3 = ℓ2 − ℓ1 ¯ d =S ℓ3 ℓ4 = d ℓ5 = ℓ3 − ℓ4 ℓ6 = 1/40 × ℓ5 ℓ7 = ℓ3 × ℓ6 ℓ8 = 0 ℓ9 = −16000 ℓ10 = if ℓ7 ≥ 0 then ℓ8 else ℓ9 ℓ11 = ℓ3 × ℓ10 ℓ12 = 20000 × ℓ3 ℓ13 = −ℓ12 − ℓ11 ℓ14 = 1 250 × ℓ13 ˙ x0 = ℓ14 ℓ15 = x0

Notation

We denote this function by fd(t, x, d) = g(t, x, d) |¯

d=Sℓ .

12 / 22

Simulink as an equation-based programming language

Kinds of equations

A Simulink model is made of four kinds of functions: the output function; (2 versions: major g and minor ˜ g) the update function of discrete-time states; the update function of continuous-time states;

ℓ1 = input ˙ x1 = ℓ15 ℓ2 = x1 ℓ3 = ℓ2 − ℓ1 ¯ d =S ℓ3 ℓ4 = d ℓ5 = ℓ3 − ℓ4 ℓ6 = 1/40 × ℓ5 ℓ7 = ℓ3 × ℓ6 ℓ8 = 0 ℓ9 = −16000 ℓ10 = if ℓ7 ≥ 0 then ℓ8 else ℓ9 ℓ11 = ℓ3 × ℓ10 ℓ12 = 20000 × ℓ3 ℓ13 = −ℓ12 − ℓ11 ℓ14 = 1 250 × ℓ13 ˙ x0 = ℓ14 ℓ15 = x0

Notation

We denote this function by fx(t, x, d) = ˜ g(t, x, d) | ˙

d=Sℓ .

12 / 22

Simulink as an equation-based programming language

Kinds of equations

A Simulink model is made of four kinds of functions: the output function; the update function of discrete-time states; the update function of continuous-time states;

Hidden equations: the 5th kind

Some blocks are associated to equations to detect zero-crossing events (only with variable step solvers).

Example

Switch block: detect when the sign of the 2nd input changes.

Notation

We denote a zero-crossing equation fz(t, x, d) .

12 / 22

Overview of numerical simulation

Goal: computing the temporal evolution of the system.

The steps of the Simulink’s simulation engine

Input: x0, d0, t0, h0; n = 0; loop until tn ≥ tend evaluate g(tn, xn, dn) ; update d′ = fd(tn, xn, dn) ; solve ˙ x(t) = fx(t, x(t), dn) over interval [tn, tn + hn] to get x(tn + hn) ; find zero crossing ; compute hn+1 ; compute tn+1 ; dn+1 = d′ ; xn+1 = x(tn + hn) ; n = n + 1 ; end loop

Remarks

Major steps: evaluation of g produces the simulation results at tn. Minor steps: evaluations of ˜ g are used as intermediate computations between tn and tn + hn.

13 / 22

Overview of numerical simulation

The temporal evolution of the system depends on a set of parameters.

Parameters (a sample)

A Simulink model is described by a set of parameters: t0= 0 start time of the simulation; tend= 10 stop time of the simulation; numerical integration methods with an absolute tolerance (atol= 10−6), a relative tolerance (rtol= 10−3); minimal (hmin= 0.2) and maximal (hmax= 5) integration step-size; zero-crossing method: adaptive or non adaptive; zero-crossing tolerance (zctol= 10 × 128 × eps).

Notation

All these parameters are represented by a record π : name → value.

Consequence

There are several semantics of Simulink.

13 / 22

Operational semantics of Simulink

· · · loop until tn ≥ tend evaluate g(tn, xn, dn) ; update d′ = fd(tn, xn, dn) ; solve ˙ x(t) = fx(t, x(t), dn) over interval [tn, tn + hn] to get x(tn + hn) ; find zero crossing ; · · · end loop

Main rules

σ(t) = π(tend) Eq, π ⊢ σ ⇒ σ simulation-end σ(t) < π(tend) Eq, π ⊢ σ

M

− → σ1 Eq, π ⊢ σ1

u

− → σ2 Eq, π ⊢ σ2

s

− → σ′ Eq, π ⊢ σ ⇒ σ′

Notation

Eq =

• g(t, x, d), fd(t, x, d), fx(t, x, d), fz(t, x, d)
• 14 / 22

Rules for output function

Recall: main rules

σ(t) < π(tend) Eq, π ⊢ σ

M

− → σ1 Eq, π ⊢ σ1

u

− → σ2 Eq, π ⊢ σ2

s

− → σ′ Eq, π ⊢ σ ⇒ σ′

Very simple rules of a classical operational semantics.

Step

Given values of continuous states x , values of discrete states d and input values. Evaluate each equations given by g(t, x, d) in the evaluation order

Sample of evaluation rules

e1, σ

• −

→ r1 e2, σ

• −

→ r2 e1 ⋄ e2, σ

• −

→ r1 ⋄ r2 arith e1, σ

• −

→ true e2, σ

• −

→ r2 if (e1, e2, e3) , σ

• −

→ r2 then

15 / 22

Rules for update discrete states

Recall: main rules

σ(t) < π(tend) Eq, π ⊢ σ

M

− → σ1 Eq, π ⊢ σ1

u

− → σ2 Eq, π ⊢ σ2

s

− → σ′ Eq, π ⊢ σ ⇒ σ′

Main idea: if the current time t is a sampling time then update the state.

• therwise keep the previous value.

Simple rules

Eq(d) = d :=S ℓ σ(t) ∈ S Eq, π ⊢ σ

u

− → σ[d → σ(ℓ)] Eq(d) = d :=S ℓ σ(t) ∈ S Eq, π ⊢ σ

u

− → σ

Remark

The numerical integration is in charge not to miss a sampling time.

16 / 22

Rules for update continuous states

Recall: Main rules

σ(t) < π(tend) Eq, π ⊢ σ

M

− → σ1 Eq, π ⊢ σ1

u

− → σ2 Eq, π ⊢ σ2

s

− → σ′ Eq, π ⊢ σ ⇒ σ′

The main features of the Simulink solver is encoded in this rules: Performing numerical integration. Handling zero-crossing. Adapting integration step-size.

Sub rules

Eq, π, σ ⊢ σ

i

− → σso Eq, π, σ ⊢ σso

zc

− → σzc Eq, π ⊢ σzc

h

− → σ′ Eq, π ⊢ σ

s

− → σ′

Remarks

Hypotheses to enforce the existence of the ODEs solution. No discrete states are updated during numerical integration. Numerical integration assumes no zero-crossing events appear.

17 / 22

Numerical integration methods – 1

There exists a huge amount of methods which may categorize in: explicit or implicit method; fixed or variable step-size; single-step or multi-step method; fixed or variable order.

Simulink’s integration methods considered

Explicit single-step fixed step-size

• de1 (Euler)
• de2 (Heun)
• de3 (Bogacki-Shampine)
• de4 (Dormand-Prince)
• de8 (Dormand-Prince)

Explicit single-step variable step-size

• de23

(Bogacki-Shampine)

• de45 (Dormand-Prince)

Remark

All these methods are members of the Runge-Kutta family which has a unified representation with Butcher Table.

18 / 22

Numerical integration methods – 2

Example

Simulink solver ode23 Mathematical description

k1 = f (tn, xn) k2 = f (tn +

1 2 hn, xn + 1 2 hk1)

k3 = f (tn +

3 4 hn, xn + 3 4 hk2)

xn+1 = xn + h

• 2

9 k1 + 1 3 k2 + 4 9 k3

• k4 = f (tn + 1 hn, xn+1 )

yn+1 = xn + h

• 7

24 k1 + 1 4 k2 + 1 3 k3 + 1 8 k4

• Butcher Table

1 2 1 2 3 4 3 4

1

2 9 1 3 4 9 2 9 1 3 4 9 7 24 1 4 1 3 1 8

Remark

The parameters set π embeds the Butcher Table of the considered integration method ⇒ handling of the several semantics of Simulink.

18 / 22

Numerical integration methods – 2

Example

Simulink solver ode23 Mathematical description

k1 = f (tn, xn) k2 = f (tn +

1 2 hn, xn + 1 2 hk1)

k3 = f (tn +

3 4 hn, xn + 3 4 hk2)

xn+1 = xn + h

• 2

9 k1 + 1 3 k2 + 4 9 k3

• k4 = f (tn + 1 hn, xn+1 )

yn+1 = xn + h

• 7

24 k1 + 1 4 k2 + 1 3 k3 + 1 8 k4

• Simulink translation

k1 = fx(tn, xn, dn); x = xn +

1 2 k1; t = tn + 1 2 hn;

k2 = fx(t, x, dn); x = xn +

3 4 k1;

· · · xn+1 = xn + h

• 2

9 k1 + 1 3 k2 + 4 9 k3

• ;

· · ·

Consequence

Our semantic rules act as an evaluation of a local expansion of the sequence of equations.

18 / 22

Numerical integration methods – 2

Evaluation rules for variable step-size solver

sc1(σM(x), Eq(x), π) , σ

• −

→ σ1 sc2(σM(x), Eq(x), π) , σ1

• −

→ σ2 check err(σ, σ1, σ2, π) = 1 Eq, π, σM ⊢ σ

i

− → σM[x → σ1(x), h → σ1(h)]

integr-success

sc1(σM(x), Eq(x), π) , σ

• −

→ σ1 sc2(σM(x), Eq(x), π) , σ1

• −

→ σ2 check err(σ, σ1, σ2, π) = 0 Eq, π, σM ⊢ σM[h → max(π(hmin), σ(h) 2 )]

i

− → σ′ Eq, π, σM ⊢ σ

i

− → σ′

integr-fail

Note: sci represents the local expansion of equations.

Solving ODEs leads to a time discretization

Equations of the form ˙ xi = ℓj are substituted by the sequence of the equations needed to realize the numerical integration.

18 / 22

Numerical integration methods – 3

Validation of the integration step

For adaptive step-size method: for all continuous state variables integration error = hn xn+1 − yn+1 ∞ max

• max
• | xn+1 |, |xn|
• , atol

rtol

• valid if

≤ rtol . Strategy: Success: go the the zero-crossing detection step. Failure: reduce the step-size hn in general only a division by 2. and restart the integration step with the new step-size.

Remark

The reduction of the step-size is done until the hmin is reached. In that case a simulation error may happen.

18 / 22

Zero-crossing event detection (non adaptive version)

Main steps

Detection of zero-crossing event Is one of the zero-crossing changed its sign between [tn, tn + hn]? Localization: if detection is true Bracket the most recent zero-crossing time using bisection method. Pass through the zero-crossing event in two steps:

Set the next major output to the left bound of the bracket time. Reset the solver with the state estimate at the right bound of bracket time.

Ingredients for the localization – 1

Linear approximation of the dynamic.

computeTz(σL, σR) = minL

• σL(t) − xL · σR(t)−σL(t)

xR−xL

: ∀x

• with minL(t) = min{ti | ti ≥ σL(t)}.

Remark: σL and σR represents the environments at the time enclosing the event.

19 / 22

Zero-crossing event detection (non adaptive version)

Main steps

Detection of zero-crossing event Is one of the zero-crossing changed its sign between [tn, tn + hn]? Localization: if detection is true Bracket the most recent zero-crossing time using bisection method. Pass through the zero-crossing event in two steps:

Set the next major output to the left bound of the bracket time. Reset the solver with the state estimate at the right bound of bracket time.

Ingredients for the localization – 2

Continuous extension (method dependent) to easily estimate state. For example, ode23

interpolX(σtn, σtn+hn, t) = (2τ 3 − 3τ 2 + 1)σtn(x) + (τ 3 − 2τ 2 + τ)(t2 − t1)p1 + (−2τ 3 + 3τ 2)σtn+hn(x) + (τ 3 − τ 2)(t2 − t1)p2

with τ = t−tn

hn

19 / 22

Update integration step-size

Goals of this step

Performance issue: increase the step-size to reduce the number of simulation loop.

• r decrease the step-size to increase the accuracy.

Scheduling issue: must not missed a discrete sampling time.

Main rule

h′ = changeH(σ(h), π) t′ = min{τ ∈ Eq(sampling) : τ > σ(t)} Eq, π, σM ⊢ σ

h

− → σ[h → min(h′, t′ − σ(t))] Note: changeH increase or decrease h in function of the integration error.

Remark

If a zero-crossing occurred the step-size is left unchanged in regards to the performance issue.

20 / 22

Conclusion

21 / 22

Conclusion and future work

Conclusion

Despite Simulink is a closed source software, the main algorithms in the solver are known.

Hope

The definition of a formal semantics will help increasing the confidence in the Simulink design flow with new verification methods.

Future work

Pursuing the formalization of the semantics: reset state, adaptive zero-crossing, integration method like ode113, etc. Development of a homemade simulator of Simulink models to validate experimentally our semantics. Development of a set-based simulator in order to address the reachability problem of Simulink models.

22 / 22