An engineering approach for the control of DPS is being developed since the eighties of the last century (Hulkó et al., 1981-2010). In the field of lumped parameters system (LPS) control, where the state/output quantities **x**(t)/**y**(t) – parameters are given as finite dimensional vectors, the actuator together with the controlled plant make up a controlled LPS. In this sense the actuators and the controlled plant as a DPS create a controlled lumped-input and distributed-parameter-output system (LDS). In this chapter the general decomposition of dynamics of controlled LDS into time and space components is introduced, which is based on numerically computed distributed parameter transient and impulse characteristics given on complex shape definition domains in 3D. Based on this decomposition a methodical framework of control synthesis decomposition into space and time tasks will be presented where in space domain approximation problems are solved and in time domain synthesis of control is realized by lumped parameter control loops. For the software support of modelling, control and design of DPS, the **Distributed Parameter Systems Blockset for MATLAB & Simulink (DPS Blockset) -** www.mathworks.com/products/connections/ has been developed within the CONNECTIONS program framework of The MathWorks, as a Third-Party Product of The MathWorks Company (Hulkó et al., 2003-2010). When solving problems in the time domain, toolboxes and blocksets of the MATLAB & Simulink software environment such as for example the Control Systems Toolbox, Simulink Control Design, System Identification Toolbox, etc. are utilized. In the space relation the approximation task is formulated as an optimization problem, where the Optimization Toolbox is made use of. A web portal named **Distributed Parameter Systems Control** www.dpscontrol.sk has been created for those interested in solving problems of DPS control (Hulkó et al., 2003-2010). This web portal features application examples from different areas of engineering practice such as the control of technological and manufacturing processes, mechatronic structures, groundwater remediation etc. Moreover this web portal offers the demo version of the **DPS Blockset** with the **Tutorial, Show, Demos** and **DPS Wizard** for download, along with the **Interactive Control** service for the interactive solution of model control problems via the

Generally in the control of lumped parameter systems the actuator and the controlled plant create a lumped parameter controlled system. In the field of DPS control the actuators together with the controlled plant - generally being a distributed-input and distributedparameter-output system (DDS) create a controlled lumped-input and distributed-

parameter-output system (LDS). Fig. 1.-3. and Fig. 6., 11., 14.

2000;… ).

Internet.

**2. DPS – DDS – LDS**

Fig. 1. Controlled DPS as LDS - heating of metal body of complex-shape

{ i }i SA (s) , { i }i SG (s) , { i }i T (x,y,z) - models of actuators DDS - distributed-input and distributed-parameter-output system { i}i U - lumped actuating quantities Ω /{ i}i Ω - complex-shape definition domain in 3D / actuation subdomains Y( x,y,z,t) – temperature field – distributed output quantity

Fig. 2. Lumped-input and distributed-parameter-output system – LDS

{ i( )}i U t – lumped input quantities

Y ,t () ( ) **x** = Y x, y,z,t – distributed output quantity

Fig. 3. General structure of lumped-input and distributed-parameter-output systems

Control of Distributed Parameter Systems - Engineering Methods

Fig. 4. i-th discrete distributed parameter impulse response of HLDS

response with maximal amplitude in point x 0,L i ∈[ ]

the means of the reduced characteristics as follows

the convolution sum *G*H x ,k U k ii i ( ) () ⊕ =

{ ii i ( )( )}q 0,k H x ,q U k q = *G* − ., see Fig. 5.

x – relation for timesteps { }k k

and Software Support in the MATLAB & Simulink Programming Environment 31

*G*H x ,k i i ( ) - partial DT impulse response in time, t - relation to the i-th input considered as

i

H (x,k) HR (x,k) H (x ,k) =

for { i i ( ) }i,k *G*H x ,k 0 ≠ , then the i-th DT distributed output quantiy in (2) can be rewritten by

At fixed i x the partial DT distributed output quantity in time direction Y x ,k is given as i i ( )

k

q 0

{ i i ( ) }q 0,k HR x ,q 1 = *G* = holds at the fixed point i x . At fixed k, the partial discrete distributed output quantity in space direction Y x,k i ( ) is given as a linear combination of elements { i ( )}q 0,k HR x,q = *G* , where the reduced discrete partial distributed characteristics

{ i ( )}q 0,k HR x,q = *G* are multiplied by corresponding elements of the set

=

i i i,k

*G G G* (4)

Y x,k H x ,k HR x,k U k i i ( ) ( ) ( ) () = ⊕ *G G* ii i (5)

ii i

H (x ,q)U (k q)

*G* − , in case the relation

{ i ( )}i,k *G*H x,k - partial DT impulse responses to the i-th input in space - x - relation

{ i ( )}i,k *G*HR x,k – reduced partial DT impulse responses to the i-th input in space,

If the reduced DT partial distributed parameter impulse responses are defined as

i

LDS - lumped-input and distributed-parameter-output system { i}i SA - actuating members of lumped input quantities { i}i GU - generators of distributed input quantities DDS - distributed-input and distributed-parameter-output system { i }i **U**(t) U (t) = - vector of lumped input quantities of LDS { i }i UA (t) - output quantities of lumped parameter actuators { i }i U ( ,t) **ξ** - distributed output quantities of generators { i}i GU U( ,t) **ξ** - overall distributed input quantity for DDS Y( ,t) Y(x, **x** = y,z,t) - distributed output quantity

Input-output dynamics of these DPS can be described, from zero initial conditions, by

$$\mathbf{Y}(\mathbf{x},\mathbf{t}) = \sum\_{i=1}^{n} \mathbf{Y}\_{i}(\mathbf{x},\mathbf{t}) = \sum\_{i=1}^{n} \mathcal{G}\_{i}(\mathbf{x},\mathbf{t}) \otimes \mathbf{U}\_{i}(\mathbf{t}) = \sum\_{i=1}^{n} \int\_{0}^{\mathbf{t}} \mathcal{G}\_{i}(\mathbf{x},\mathbf{t}) \mathbf{U}\_{i}(\mathbf{t} - \mathbf{\tau}) \,\mathrm{d}\mathbf{\tau} \tag{1}$$

or in discrete form

$$\mathbf{Y}(\mathbf{x},\mathbf{k}) = \sum\_{i=1}^{n} \mathbf{Y}\_{i}(\mathbf{x},\mathbf{k}) = \sum\_{i=1}^{n} \mathbf{G} \mathbf{H}\_{i}(\mathbf{x},\mathbf{k}) \oplus \mathbf{U}\_{i}(\mathbf{k}) = \sum\_{i=1}^{n} \sum\_{\mathbf{q}=0}^{k} \mathbf{G} \mathbf{H}\_{i}(\mathbf{x},\mathbf{q}) \mathbf{U}\_{i}(\mathbf{k}-\mathbf{q}) \tag{2}$$

where ⊗ marks convolution product and ⊕ marks convolution sum, *G*i(**x**,t) – distributed parameter impulse response of LDS to the i-th input, *G*Hi(**x**,k) – discrete time (DT) distributed parameter impulse response of LDS with zero-order hold units H - HLDS to the i-th input, Yi(**x**,t) - distributed output quantity of LDS to the i-th input, Yi(**x**,k) – DT distributed output quantity of HLDS to the i-th input. For simplicity in this chapter distributed quantities are considered mostly as continuous scalar quantity fields with unit sampling interval in time domain. Whereas DT distributed parameter step responses { i( )}i,k *H*H ,k **x** of HLDS can be computed by common analytical or numerical methods then DT distributed parameter impulse responses can be obtained as

$$\left\{\mathcal{G}\mathbb{H}\_{i}\left(\mathbf{x},\mathbf{k}\right) = \mathcal{H}\mathbb{H}\_{i}\left(\mathbf{x},\mathbf{k}\right) - \mathcal{H}\mathbb{H}\_{i}\left(\mathbf{x},\mathbf{k}\cdot\mathbf{1}\right)\right\}\_{i,\mathbf{k}}\tag{3}$$

#### **3. Decomposition of dynamics**

The process of dynamics decomposition shall be started from DT distributed parameter step and impulse responses of the analysed LDS. For an illustration, procedure of decomposition of dynamics and control synthesis will be shown on the LDS with zero-order hold units H – HLDS – distributed only on the interval [0,L , with output quantity ]

() () n i i 1 Y x,k Y x,k = = discretised in time relation and continuous in space relation on this

interval. Nevertheless the following results are valid in general both for continous or discrete distributed quantities in space relation given on compex-shape definition domains over 3D as well.

Input-output dynamics of these DPS can be described, from zero initial conditions, by

i 1 i 1 i 1

n n n k

i 1 i 1 i 1q 0

= = = =

DT distributed parameter impulse responses can be obtained as

**3. Decomposition of dynamics**

() () n i i 1 Y x,k Y x,k =

over 3D as well.

== =

( )( ) nn n t i ii i i 0

i i i ii

Y( ,k) Y ( ,k) H ( ,k) U (k) H ( ,q)U (k q)

where ⊗ marks convolution product and ⊕ marks convolution sum, *G*i(**x**,t) – distributed parameter impulse response of LDS to the i-th input, *G*Hi(**x**,k) – discrete time (DT) distributed parameter impulse response of LDS with zero-order hold units H - HLDS to the i-th input, Yi(**x**,t) - distributed output quantity of LDS to the i-th input, Yi(**x**,k) – DT distributed output quantity of HLDS to the i-th input. For simplicity in this chapter distributed quantities are considered mostly as continuous scalar quantity fields with unit sampling interval in time domain. Whereas DT distributed parameter step responses { i( )}i,k *H*H ,k **x** of HLDS can be computed by common analytical or numerical methods then

The process of dynamics decomposition shall be started from DT distributed parameter step and impulse responses of the analysed LDS. For an illustration, procedure of decomposition of dynamics and control synthesis will be shown on the LDS with zero-order hold units H – HLDS – distributed only on the interval [0,L , with output quantity ]

interval. Nevertheless the following results are valid in general both for continous or discrete distributed quantities in space relation given on compex-shape definition domains

= discretised in time relation and continuous in space relation on this

= = ⊗ = τ −τ τ **xxx x** *G G* (1)

**xx x** = = ⊕= *G G* **x** − (2)

{ iii () () ( )}i,k *GHH* H ,k H ,k H ,k-1 **xxx** = − (3)

Y( ,t) Y ( ,t) ( ,t) U (t) , U t d

LDS - lumped-input and distributed-parameter-output system { i}i SA - actuating members of lumped input quantities { i}i GU - generators of distributed input quantities

DDS - distributed-input and distributed-parameter-output system

{ i }i **U**(t) U (t) = - vector of lumped input quantities of LDS { i }i UA (t) - output quantities of lumped parameter actuators { i }i U ( ,t) **ξ** - distributed output quantities of generators { i}i GU

U( ,t) **ξ** - overall distributed input quantity for DDS Y( ,t) Y(x, **x** = y,z,t) - distributed output quantity

or in discrete form

Fig. 4. i-th discrete distributed parameter impulse response of HLDS *G*H x ,k i i ( ) - partial DT impulse response in time, t - relation to the i-th input considered as response with maximal amplitude in point x 0,L i ∈[ ] { i ( )}i,k *G*H x,k - partial DT impulse responses to the i-th input in space - x - relation { i ( )}i,k *G*HR x,k – reduced partial DT impulse responses to the i-th input in space, x – relation for timesteps { }k k If the reduced DT partial distributed parameter impulse responses are defined as

$$\left\{\mathcal{G}\text{H}\mathbf{R}\_{\text{i}}(\mathbf{x},\mathbf{k})=\frac{\mathcal{G}\mathbf{H}\_{\text{i}}(\mathbf{x},\mathbf{k})}{\mathcal{G}\mathbf{H}\_{\text{i}}(\mathbf{x}\_{\text{i}},\mathbf{k})}\right\}\_{\text{i},\mathbf{k}}\tag{4}$$

for { i i ( ) }i,k *G*H x ,k 0 ≠ , then the i-th DT distributed output quantiy in (2) can be rewritten by the means of the reduced characteristics as follows

$$\mathbf{Y}\_{i}(\mathbf{x},\mathbf{k}) = \mathfrak{G}\mathbf{H}\_{i}(\mathbf{x}\_{i},\mathbf{k})\mathfrak{G}\mathbf{H}\mathbf{R}\_{i}(\mathbf{x},\mathbf{k})\oplus\mathbf{U}\_{i}(\mathbf{k})\tag{5}$$

At fixed i x the partial DT distributed output quantity in time direction Y x ,k is given as i i ( ) the convolution sum *G*H x ,k U k ii i ( ) () ⊕ = k ii i q 0 H (x ,q)U (k q) = *G* − , in case the relation { i i ( ) }q 0,k HR x ,q 1 = *G* = holds at the fixed point i x . At fixed k, the partial discrete distributed output quantity in space direction Y x,k i ( ) is given as a linear combination of elements { i ( )}q 0,k HR x,q = *G* , where the reduced discrete partial distributed characteristics { i ( )}q 0,k HR x,q = *G* are multiplied by corresponding elements of the set { ii i ( )( )}q 0,k H x ,q U k q = *G* − ., see Fig. 5.

This decomposition is valid for all given lumped input { i}i U and corresponding output quantities { i ( )}i,k Y x,k - thus we obtain time and space components of HLDS dynamics in the following form:

*Time Components of Dynamics* { i i ( )}i,k *G*H x ,k – for given i and chosen xi - variable k *Space Components of Dynamics* { i( )}i,k *G*HR x,k – for given i and chosen k – variable x Also reduced components of single distributed output quantities are

$$\left\{ \mathbf{Y} \mathbf{R}\_{\mathbf{i}} \left( \mathbf{x}, \mathbf{k} \right) \right\}\_{\mathbf{i}} = \left\{ \frac{\mathbf{Y}\_{\mathbf{i}} \left( \mathbf{x}, \mathbf{k} \right)}{\mathbf{Y}\_{\mathbf{i}} \left( \mathbf{x}\_{\mathbf{i}}, \mathbf{k} \right)} \right\}\_{\mathbf{i}, \mathbf{k}} \tag{6}$$

Control of Distributed Parameter Systems - Engineering Methods

**4. Distributed parameter systems of control**

For the steady-state, when k → ∞

then in the steady-state

space relation as well.

**4.1 Open-loop control**

E x, ( ) ∞ in the quadratic norm:

lumped parameter systems of control.

and Software Support in the MATLAB & Simulink Programming Environment 33

() ( ) () n n ii i ii i

i 1 i 1

= =

Y(x, ) Y (x , )YR x, Y x , HR x,

When distributed quantities are used in discrete form as finite sequences of quantities, the discretization in space domain is usually considered by the computational nodes of the numerical model of the controlled DPS over the compex-shape definition domain in 3D.

Based on decomposition of HLDS dynamics into time and space components, possibilities for control synthesis are also suggested by an analogous approach. In this section a methodical framework for the decomposition of control synthesis into space and time problems will be presented by select demonstration control problems. In the space domain control synthesis will be solved as a sequence of approximation tasks on the set of space components of controlled system dynamics, where distributed parameter quantities in particular sampling times are considered as continuous functions on the interval [0,L as ] elements of strictly convex normed linear space **X** with quadratic norm. It is necessary to note as above that the following results are valid in general for DPS given on compexshape definition domains in 3D both for continous or discrete distributed quantities, in the

In the time domain, the control synthesis solutions are based on synthesis methods of DT

Assume the open-loop control of a distributed parameter system, where dynamic characteristics give an ideal representation of controlled system dynamics and V x,t 0 ( ) = , that is with zero initial steady-state, in which all variables involved are equal to zero – see see Fig. 6 for reference. Let us consider a step change of distributed reference quantity - W x,k W x, () ( ) = ∞ , see Fig. 7. For simplicity let the goal of the control synthesis is to generate a sequence of control inputs U k( ) in such fashion that in the steady-state, for k → ∞ , the control error E x,k W x, Y x,k () ( )() = ∞− will approach its minimal value

min E x, min W x, Y x, E x, ( ) ∞ = ∞− ∞ = ∞ ()() ( )

() ( ) () n ii i

∞− ∞ ∞ *H* (10)

First, an approximation problem will be solved in the space synthesis (SS) block:

i 1 min W x, W x , HR x, =

(9)

{ i i () ()}i,k YR x, HR x, ∞→ ∞ *H* (7)

∞= ∞ ∞= ∞ ∞ *H* (8)

then { i i ( ) }i,k Y x ,k 0 ≠ can be considered as time components and { i( )}i,k YR x,k as space components of the output quantities.

When reduced steady-state distributed parameter transient responses are introduced { i ii ( )} { () ( ) i } i i *H HH* HR x, H x, / H x , ∞= ∞ ∞ - for { i i ( ) }i *H*H x, 0 ∞ ≠ - and discrete transfer functions { i i ( )}i SH x ,z are assigned to partial distributed parameter transient responses with maximal amplitudes at points { i}i x on the interval [0,L , we obtain time and space ] components of HLDS dynamics for steady-state as:

*Time Components of Dynamics* { i i ( )}i SH x ,z - for given i and chosen xi - variable z *Space Components of Dynamics* { i ( )}i *H*HR x,∞ - for given i in ∞ – variable x

Fig. 5. Partial distributed output quantities in time and space direction

Ui - i-th DT lumped input quantity

Y x ,k - i-th partial DT distributed output qu i i ( ) antity in time domain at chosen point i x

Y x,k i ( ) / YR x,k i ( ) - i-th partial distributed output/reduced output quantity in space direction

For the steady-state, when k → ∞

$$\left\{ \mathbf{Y} \mathsf{R}\_{\mathrm{i}} (\mathsf{x}\_{\mathsf{'}} \mathsf{s}\_{\mathsf{'}}) \to \mathsf{M} \mathsf{H} \mathsf{R}\_{\mathrm{i}} (\mathsf{x}\_{\mathsf{'}} \mathsf{s}\_{\mathsf{'}}) \right\}\_{\mathsf{i},\mathsf{k}} \tag{7}$$

then in the steady-state

32 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

This decomposition is valid for all given lumped input { i}i U and corresponding output quantities { i ( )}i,k Y x,k - thus we obtain time and space components of HLDS dynamics in

*Time Components of Dynamics* { i i ( )}i,k *G*H x ,k – for given i and chosen xi - variable k

*Space Components of Dynamics* { i( )}i,k *G*HR x,k – for given i and chosen k – variable x

i i

*Time Components of Dynamics* { i i ( )}i SH x ,z - for given i and chosen xi - variable z

*Space Components of Dynamics* { i ( )}i *H*HR x,∞ - for given i in ∞ – variable x

Fig. 5. Partial distributed output quantities in time and space direction

Y x ,k - i-th partial DT distributed output qu i i ( ) antity in time domain at chosen point i x

Y x,k i ( ) / YR x,k i ( ) - i-th partial distributed output/reduced output quantity in space

Ui - i-th DT lumped input quantity

direction

{ } ( ) ( )

Y x,k YR x,k Y x ,k

then { i i ( ) }i,k Y x ,k 0 ≠ can be considered as time components and { i( )}i,k YR x,k as space

When reduced steady-state distributed parameter transient responses are introduced { i ii ( )} { () ( ) i } i i *H HH* HR x, H x, / H x , ∞= ∞ ∞ - for { i i ( ) }i *H*H x, 0 ∞ ≠ - and discrete transfer functions { i i ( )}i SH x ,z are assigned to partial distributed parameter transient responses with maximal amplitudes at points { i}i x on the interval [0,L , we obtain time and space ]

( ) i

i i i,k

= (6)

Also reduced components of single distributed output quantities are

the following form:

components of the output quantities.

components of HLDS dynamics for steady-state as:

$$\mathbf{Y}(\mathbf{x},\boldsymbol{\omega}\circ) = \sum\_{i=1}^{n} \mathbf{Y}\_{i}(\mathbf{x}\_{i},\boldsymbol{\omega}\circ) \mathbf{Y} \mathbf{R}\_{i}\left(\mathbf{x},\boldsymbol{\omega}\circ\right) = \sum\_{i=1}^{n} \mathbf{Y}\_{i}\left(\mathbf{x}\_{i},\boldsymbol{\omega}\circ\right) \mathcal{H} \mathbf{H} \mathbf{R}\_{i}\left(\mathbf{x},\boldsymbol{\omega}\circ\right) \tag{8}$$

When distributed quantities are used in discrete form as finite sequences of quantities, the discretization in space domain is usually considered by the computational nodes of the numerical model of the controlled DPS over the compex-shape definition domain in 3D.

#### **4. Distributed parameter systems of control**

Based on decomposition of HLDS dynamics into time and space components, possibilities for control synthesis are also suggested by an analogous approach. In this section a methodical framework for the decomposition of control synthesis into space and time problems will be presented by select demonstration control problems. In the space domain control synthesis will be solved as a sequence of approximation tasks on the set of space components of controlled system dynamics, where distributed parameter quantities in particular sampling times are considered as continuous functions on the interval [0,L as ] elements of strictly convex normed linear space **X** with quadratic norm. It is necessary to note as above that the following results are valid in general for DPS given on compexshape definition domains in 3D both for continous or discrete distributed quantities, in the space relation as well.

In the time domain, the control synthesis solutions are based on synthesis methods of DT lumped parameter systems of control.

#### **4.1 Open-loop control**

Assume the open-loop control of a distributed parameter system, where dynamic characteristics give an ideal representation of controlled system dynamics and V x,t 0 ( ) = , that is with zero initial steady-state, in which all variables involved are equal to zero – see see Fig. 6 for reference. Let us consider a step change of distributed reference quantity - W x,k W x, () ( ) = ∞ , see Fig. 7. For simplicity let the goal of the control synthesis is to generate a sequence of control inputs U k( ) in such fashion that in the steady-state, for k → ∞ , the control error E x,k W x, Y x,k () ( )() = ∞− will approach its minimal value E x, ( ) ∞ in the quadratic norm:

$$\min \left\| \mathbb{E} \left( \mathbf{x}\_{\prime \prime \ast} \right) \right\| = \min \left\| \mathbb{W} \left( \mathbf{x}\_{\prime \prime \ast} \right) - \mathbb{Y} \left( \mathbf{x}\_{\prime \prime \ast} \right) \right\| = \left\| \bar{\mathbb{E}} \left( \mathbf{x}\_{\prime \prime \ast} \right) \right\| \tag{9}$$

First, an approximation problem will be solved in the space synthesis (SS) block:

$$\min \left\| \mathbf{W} \left( \mathbf{x}\_{\prime \prime \ast} \right) - \sum\_{i=1}^{n} \mathbf{W}\_{i} \left( \mathbf{x}\_{i \prime \ast} \right) \mathcal{H} \mathbf{H} \mathbf{R}\_{i} \left( \mathbf{x}\_{\prime \prime \ast} \right) \right\| \tag{10}$$

Control of Distributed Parameter Systems - Engineering Methods

n

min W x, W x , HR x,

i 1 n

=

i 1

=

Fig. 7. Step change of distributed reference quantity HLDS - controlled system with zero-order hold units

Fig. 8. Solution of the approximation problem

{ i}i U - lumped control quantities

HLDS - controlled system with zero-order hold units

{ i ii } { ( )} i i W W x, = ∞ - optimal approximation parameters, lumped references

{ i ( )}i *H*HR x,∞ - reduced steady-state distributed parameter transient responses

W x,k W () ( ) = ∞ x, - step change of distributed reference quantity

{ i}i U - lumped control quantities

and Software Support in the MATLAB & Simulink Programming Environment 35

() ( ) () () ()

(11)

W x, W x , HR x, W x, WO x,

= ∞− ∞ ∞ = ∞− ∞

() ( ) ()

ii i

ii i

*H*

∞− ∞ ∞ =

*H*

Fig. 6. Distributed parameter open-loop system of control

LDS - lumped-input and distributed-parameter-output system

H - zero-order hold units

HLDS - controlled system with zero-order hold units

CS - control synthesis

TS - time part of control synthesis

SS - space part of control synthesis

Y x,t /W x,k W x, () ( ) ( ) = ∞ - distributed controlled/reference quantity

V x,t ( ) - distributed disturbance quantity

W ={ i i ( )}i W x,∞ - vector of lumped reference quantities

U k( ) - vector of lumped control quantities

{ ii i ( ) ()}i SH x ,z / HR x, *H* ∞ - time/space components of controlled system dynamics

where { i ( )}i *H*HR x,∞ are steady-state reduced distributed parameter transient responses of the controlled system – HLDS and { i i ( )} W x, i ∞ are parameters of approximation. Functions

{ i ( )}i *H*HR x,∞ form a finite-dimensional subspace of approximation functions *Fn* in the strictly convex normed linear space of distributed parameter quantities **X** on [0,L with ] quadratic norm, where the approximation problem is to be solved, see Fig. 8. for reference. From approximation theory in this relation is known the theorem:

*Let Fn be a finite-dimensional subspace of a strictly convex normed linear space X. Then, for each f*∈ *X, there exists a unique element of best approximation.*

(Shadrin, 2005). So the solution of the approximation problem (10) is guaranteed as a unique

best approximation () ( ) () n ii i i 1 WO x, W x , HR x, = ∞= ∞ ∞ *H* , where { i i ( )}i W W x, = ∞ is the

vector of optimal approximation parameters. Hence:

Fig. 6. Distributed parameter open-loop system of control

HLDS - controlled system with zero-order hold units

H - zero-order hold units

TS - time part of control synthesis SS - space part of control synthesis

V x,t ( ) - distributed disturbance quantity

U k( ) - vector of lumped control quantities

CS - control synthesis

W

LDS - lumped-input and distributed-parameter-output system

Y x,t /W x,k W x, () ( ) ( ) = ∞ - distributed controlled/reference quantity

{ ii i ( ) ()}i SH x ,z / HR x, *H* ∞ - time/space components of controlled system dynamics

where { i ( )}i *H*HR x,∞ are steady-state reduced distributed parameter transient responses of the controlled system – HLDS and { i i ( )} W x, i ∞ are parameters of approximation. Functions { i ( )}i *H*HR x,∞ form a finite-dimensional subspace of approximation functions *Fn* in the strictly convex normed linear space of distributed parameter quantities **X** on [0,L with ] quadratic norm, where the approximation problem is to be solved, see Fig. 8. for reference.

*Let Fn be a finite-dimensional subspace of a strictly convex normed linear space X. Then, for each*

(Shadrin, 2005). So the solution of the approximation problem (10) is guaranteed as a unique

∞= ∞ ∞ *H* , where { i i ( )}i W W x, = ∞

is the

ii i

={ i i ( )}i W x,∞ - vector of lumped reference quantities

From approximation theory in this relation is known the theorem:

i 1 WO x, W x , HR x, =

*f*∈ *X, there exists a unique element of best approximation.*

vector of optimal approximation parameters. Hence:

best approximation () ( ) () n

$$\begin{aligned} &\min \left\| \mathbf{W} \left( \mathbf{x}\_{\prime \prime \circ} \right) - \sum\_{i=1}^{n} \mathbf{W}\_{i} \left( \mathbf{x}\_{i} \, \circ \right) \mathcal{H} \mathbf{R}\_{i} \left( \mathbf{x}\_{\prime \prime \circ} \right) \right\| = \\ &= \left\| \mathbf{W} \left( \mathbf{x}\_{\prime \prime \circ} \right) - \sum\_{i=1}^{n} \bar{\mathbf{W}}\_{i} \left( \mathbf{x}\_{i} \, \circ \right) \mathcal{H} \mathbf{R}\_{i} \left( \mathbf{x}\_{\prime \prime \circ} \right) \right\| = \left\| \mathbf{W} \left( \mathbf{x}\_{\prime \prime \circ} \right) - \bar{\mathbf{W}} \mathbf{O} \left( \mathbf{x}\_{\prime \prime \circ} \right) \right\| \end{aligned} \tag{11}$$

Fig. 7. Step change of distributed reference quantity

HLDS - controlled system with zero-order hold units { i}i U - lumped control quantities

W x,k W () ( ) = ∞ x, - step change of distributed reference quantity

Fig. 8. Solution of the approximation problem

HLDS - controlled system with zero-order hold units { i}i U - lumped control quantities { i ii } { ( )} i i W W x, = ∞ - optimal approximation parameters, lumped references { i ( )}i *H*HR x,∞ - reduced steady-state distributed parameter transient responses

Control of Distributed Parameter Systems - Engineering Methods

{ i ( )}i,k U k - lumped control quantities

for k → ∞ converging to reference quantities

see Fig. 9. – 10. for reference.

{ i i ( )}i,k E x ,k

W k /U k ( ) ( )

{ i i ( )} { i ii } { ( )} i,k i i Y x ,k / W W x , = ∞ - controlled/reference quantities

When the individual components of the vector { i i ( )}i W W x, = ∞

and Software Support in the MATLAB & Simulink Programming Environment 37

control loops: { ii i ( ) ()}i SH x ,z ,R z , the control processes take place. The applied control laws result in the sequences of control inputs: { i ( )}i,k U k , and respectively the output quantities,

Values of these lumped controlled quantities in new steady-state will be further denoted as

Fig. 10. Quantities of distributed parameter open-loop control in new steady-state

{ ii i ii () ()}i Y x , /W W x , ∞= ∞ - controlled/reference quantities in new steady-state

HLDS - controlled system with zero-order hold units

Y x, ( ) ∞ - controlled distributed quantity in new steady-state

{ i}i U - lumped control quantities

{ i i ( )}i E x,∞ - lumped control errors

are input in the particular

{ ii ii ii () ( ) ( )}i,k Y x ,k Y x , W x , → ∞= ∞ (13)

{ ii ii () ()}i Y x, W x, ∞= ∞ (14)

W x, ( ) ∞ - distributed reference quantity

WO x, ( ) ∞ - unique best approximation of reference quantity

Let us assume vector W enters the block of time synthesis (TS). In this block, there are "n" single-input /single-output (SISO) lumped parameter control loops: { ii i ( ) ()}i SH x ,z ,R z , see Fig. 9. for reference. The controlled systems of these loops are lumped parameter systems assigned to HLDS as time components of dynamics: { i i ( )}i SH x ,z . Controllers, { i ( )}i R z , are to be chosen such that for k → ∞ the following relation holds:

$$\mathbb{E}\left[\lim\_{\mathbf{k}\to\boldsymbol{\ast}}\breve{\mathbf{E}}\_{\mathbf{i}}\left(\mathbf{x}\_{\mathbf{i}},\mathbf{k}\right)=\lim\_{\mathbf{k}\to\boldsymbol{\ast}}\bigg[\breve{W}\_{\mathbf{i}}\left(\mathbf{x}\_{\mathbf{i}},\boldsymbol{\ast}\right)-\varUpsilon\_{\mathbf{i}}\left(\mathbf{x}\_{\mathbf{i}},\mathbf{k}\right)\bigg]=0\bigg]\_{\mathbf{i},\mathbf{k}}\tag{12}$$

Fig. 9. SISO lumped parameter control loops in the block TS

TS - time part of control synthesis

{ i i ( )}i SH x ,z - time components of HLDS dynamics

{ i ( )}i R z - lumped parameter controllers

{ i i ( )} { i ii } { ( )} i,k i i Y x ,k / W W x , = ∞ - controlled/reference quantities { i i ( )}i,k E x ,k - lumped control errors { i ( )}i,k U k - lumped control quantities

W k /U k ( ) ( ) - vector of lumped reference/control quantities

When the individual components of the vector { i i ( )}i W W x, = ∞ are input in the particular control loops: { ii i ( ) ()}i SH x ,z ,R z , the control processes take place. The applied control laws result in the sequences of control inputs: { i ( )}i,k U k , and respectively the output quantities, for k → ∞ converging to reference quantities

$$\left\{ \mathbf{Y}\_{\mathrm{i}} \left( \mathbf{x}\_{\mathrm{i}}, \mathbf{k} \right) \to \mathbf{Y}\_{\mathrm{i}} \left( \mathbf{x}\_{\mathrm{i}}, \rightsquigarrow \right) = \bar{\mathbf{W}}\_{\mathrm{i}} \left( \mathbf{x}\_{\mathrm{i}}, \rightsquigarrow \right) \right\}\_{\mathrm{i}, \mathrm{k}} \tag{13}$$

Values of these lumped controlled quantities in new steady-state will be further denoted as

$$\left\{ \bar{\mathbf{Y}}\_{\mathbf{i}} \left( \mathbf{x}\_{\mathbf{i}}, \leadsto \right) = \bar{\mathbf{W}}\_{\mathbf{i}} \left( \mathbf{x}\_{\mathbf{i}}, \leadsto \right) \right\}\_{\mathbf{i}} \tag{14}$$

see Fig. 9. – 10. for reference.

36 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

single-input /single-output (SISO) lumped parameter control loops: { ii i ( ) ()}i SH x ,z ,R z , see Fig. 9. for reference. The controlled systems of these loops are lumped parameter systems assigned to HLDS as time components of dynamics: { i i ( )}i SH x ,z . Controllers, { i ( )}i R z , are

> { i i () ( )() ii ii } k k i,k lim E x ,k lim W x , Y x ,k 0

= ∞− =

enters the block of time synthesis (TS). In this block, there are "n"

(12)

W x, ( ) ∞ - distributed reference quantity

Let us assume vector W

WO x, ( ) ∞ - unique best approximation of reference quantity

to be chosen such that for k → ∞ the following relation holds:

Fig. 9. SISO lumped parameter control loops in the block TS

{ i i ( )}i SH x ,z - time components of HLDS dynamics

TS - time part of control synthesis

{ i ( )}i R z - lumped parameter controllers

→∞ →∞

{ i}i U - lumped control quantities { ii i ii () ()}i Y x , /W W x , ∞= ∞ - controlled/reference quantities in new steady-state { i i ( )}i E x,∞ - lumped control errors Y x, ( ) ∞ - controlled distributed quantity in new steady-state

Control of Distributed Parameter Systems - Engineering Methods

time components of partial output quantities { i i ( )}i Y x ,k

{ ii ii ii () ( )()}i E x ,k W x , Y x ,k = ∞− . These sequences { i i ( )}i E x ,k

and give { i ( )}i U k , which enter into HLDS with { i i ( )}i Y x ,k

parameter control loops { ii i ( ) ()}i SH x ,z ,R z , Fig. 9.

subspace of { i ( )}i,k YR x,k

U k and { i i ( )}i Y x ,k

components of output quantity

{ i ( )}i

and Software Support in the MATLAB & Simulink Programming Environment 39

generated. Block TS in Fig. 12., contains the controllers{ i ( )}i R z designed as the controllers for SISO lumped parameter control loops { ii i ( ) ()}i SH x ,z ,R z with respect to the request formulated by equation (12). In the k-th step in block SS2 at approximation of Y x,k ( ) on the

> () ( ) () n ii i

− (19)

there are relations { i i ( )}i SH x ,z . - This analysis of control

are obtained, in block SS1

enter into the TS on{ i ( )}i

on the SS2 output - among

R z

i 1 min Y x,k Y x ,k YR x,k =

reference quantities { i i ( )}i W x,∞ are computed. Then on the output of the algebraic block is

synthesis process shows that synthesis in time domain is realized on the level of one

Fig. 11. Distributed parameter closed-loop system of control with reduced space

W x, ( ) ∞ - distributed reference quantity

WO x, ( ) ∞ - unique best approximation of reference quantity

E x, ( ) ∞ - distributed control error with minimal norm

Then according to equations (12-14) for the new steady-state it holds: { ii i ii i ( ) () ( ) ()}i Y x , HR x, W x , HR x, ∞ ∞= ∞ ∞ *H H* , which implies that the overall distributed output quantity at the time k → ∞ : Y x, ( ) ∞ gives the unique best approximation of the distributed reference variable: W x, ( ) ∞

$$\bar{\mathbf{Y}}(\mathbf{x},\boldsymbol{\omega}\circ\rangle = \sum\_{i=1}^{n} \bar{\mathbf{Y}}\_{i}(\mathbf{x}\_{i},\boldsymbol{\omega}\circ\rangle \mathcal{H}\mathbf{H}\mathbf{R}\_{i}(\mathbf{x},\boldsymbol{\omega}\circ\rangle = \sum\_{i=1}^{n} \bar{\mathbf{W}}\_{i}(\mathbf{x}\_{i},\boldsymbol{\omega}\circ\rangle \mathcal{H}\mathbf{H}\mathbf{R}\_{i}(\mathbf{x},\boldsymbol{\omega}\circ\rangle = \bar{\mathbf{W}}\mathbf{O}\left(\mathbf{x},\boldsymbol{\omega}\circ\right) \tag{15}$$

Therefore the control error has a unique form as well, with minimal quadratic norm

$$\begin{aligned} \left\| \bar{\mathcal{E}} \left( \mathbf{x}, \curvearrowright \right) \right\| &= \left\| \mathcal{W} \left( \mathbf{x}, \curvearrowright \right) - \bar{\mathcal{Y}} \left( \mathbf{x}, \curvearrowright \right) \right\| = \left\| \mathcal{W} \left( \mathbf{x}, \curvearrowright \right) - \sum\_{i=1}^{n} \bar{\mathcal{Y}}\_{i} \left( \mathbf{x}\_{i}, \curvearrowright \right) \mathcal{H} \mathrm{IR}\_{i} \left( \mathbf{x}, \curvearrowright \right) \right\| \\ &= \left\| \mathcal{W} \left( \mathbf{x}, \curvearrowright \right) - \sum\_{i=1}^{n} \bar{\mathcal{W}}\_{i} \left( \mathbf{x}\_{i}, \curvearrowright \right) \mathcal{H} \mathrm{IR}\_{i} \left( \mathbf{x}, \curvearrowright \right) \right\| = \left\| \mathcal{W} \left( \mathbf{x}, \curvearrowright \right) - \bar{\mathcal{W}} \mathcal{O} \left( \mathbf{x}, \curvearrowright \right) \right\| \end{aligned} \tag{16}$$

Thus the control task, defined at equation (9), is accomplished with { i i ( ) }i E x, 0 ∞ = - see Fig. 10. for reference. As the conclusion of this section we may state that the control synthesis in open-loop control system is realized as:

*Time Tasks of Control Synthesis* – in lumped parameter control loops

*Space Tasks of Control Synthesis* – as approximation task.

When mathematical models cannot provide an ideal representation of controlled DPS dynamics and disturbances are present with an overall effect on the output in steady-state expressed by EY x, ( ) ∞ then

$$\begin{aligned} \left\|\mathbf{W}(\mathbf{x},\boldsymbol{\omega}\circ)-\bar{\mathbf{Y}}(\mathbf{x},\boldsymbol{\omega}\circ)\pm\mathrm{EY}(\mathbf{x},\boldsymbol{\omega}\circ)\right\| &\leq \left\|\mathbf{W}(\mathbf{x},\boldsymbol{\omega}\circ)-\bar{\mathbf{Y}}(\mathbf{x},\boldsymbol{\omega}\circ)\right\| + \left\|\mathrm{EY}(\mathbf{x},\boldsymbol{\omega}\circ)\right\| = \\ \left\|\mathbf{W}(\mathbf{x},\boldsymbol{\omega}\circ)-\bar{\mathbf{W}}\mathbf{O}(\mathbf{x},\boldsymbol{\omega}\circ)\right\| + \left\|\mathrm{EY}(\mathbf{x},\boldsymbol{\omega}\circ)\right\| &= \left\|\bar{\mathbf{E}}(\mathbf{x},\boldsymbol{\omega}\circ)\right\| + \left\|\mathrm{EY}(\mathbf{x},\boldsymbol{\omega}\circ)\right\| \end{aligned} \tag{17}$$

Finally at the design stage of a control system, for a given desired quality of control δ in space domain, it is necessary to choose appropriate number and layout of actuators for the fulfillment of this requirement

$$\left\|\bar{\mathbf{E}}\left(\mathbf{x}\_{\prime}\approx\right)\right\| + \left\|\mathbf{E}Y(\mathbf{x}\_{\prime}\approx)\right\| \leq \delta \tag{18}$$

#### **4.2 Closed-loop control with block RHLDS**

Let us consider now a distributed parameter feedback control loop with initial conditions identical as the case above, see Fig. 11. In blocks SS1 a SS2 approximation problems are solved while in block RHLDS reduced distributed output quantities { i( )}i,k YR x,k are

Then according to equations (12-14) for the new steady-state it holds: { ii i ii i ( ) () ( ) ()}i Y x , HR x, W x , HR x, ∞ ∞= ∞ ∞ *H H* , which implies that the overall distributed output quantity at the time k → ∞ : Y x, ( ) ∞ gives the unique best approximation of the

> () ( ) () ( ) () () n n ii i ii i

( ) ()() () ( ) ()

Thus the control task, defined at equation (9), is accomplished with { i i ( ) }i E x, 0 ∞ = - see Fig. 10. for reference. As the conclusion of this section we may state that the control synthesis in

When mathematical models cannot provide an ideal representation of controlled DPS dynamics and disturbances are present with an overall effect on the output in steady-state

()() () ()() ( )

Finally at the design stage of a control system, for a given desired quality of control δ in space domain, it is necessary to choose appropriate number and layout of actuators for the

Let us consider now a distributed parameter feedback control loop with initial conditions identical as the case above, see Fig. 11. In blocks SS1 a SS2 approximation problems are solved while in block RHLDS reduced distributed output quantities { i( )}i,k YR x,k are

∞− ∞± ∞ ≤ ∞− ∞ + ∞ =

(17)

E x, EY x, ( ) ∞ + ∞ ≤δ ( ) (18)

() () ( ) ( ) ( )

∞− ∞ + ∞ = ∞ + ∞

W x, WO x, EY x, E x, EY x,

W x, Y x, EY x, W x, Y x, EY x,

∞ = ∞− ∞ = ∞− ∞ ∞ =

E x, W x, Y x, W x, Y x , HR x,

∞ ∞ ∞= ∞ ∞= ∞ = *H H* (15)

n

i 1

=

ii i

*H*

(16)

Y x, Y x , HR x, W x , HR x, WO x,

Therefore the control error has a unique form as well, with minimal quadratic norm

() ( ) () () ()

∞− ∞ ∞ = ∞− ∞

W x, W x , HR x, W x, WO x,

W x, ( ) ∞ - distributed reference quantity

distributed reference variable: W x, ( ) ∞

n

i 1

=

open-loop control system is realized as:

expressed by EY x, ( ) ∞ then

fulfillment of this requirement

**4.2 Closed-loop control with block RHLDS**

WO x, ( ) ∞ - unique best approximation of reference quantity

i 1 i 1

= =

ii i

*Time Tasks of Control Synthesis* – in lumped parameter control loops

*Space Tasks of Control Synthesis* – as approximation task.

*H*

E x, ( ) ∞ - distributed control error with minimal norm

generated. Block TS in Fig. 12., contains the controllers{ i ( )}i R z designed as the controllers for SISO lumped parameter control loops { ii i ( ) ()}i SH x ,z ,R z with respect to the request formulated by equation (12). In the k-th step in block SS2 at approximation of Y x,k ( ) on the subspace of { i ( )}i,k YR x,k

$$\min \left\| \mathbf{Y}(\mathbf{x}, \mathbf{k}) - \sum\_{i=1}^{n} \mathbf{Y}\_{i}(\mathbf{x}\_{i}, \mathbf{k}) \mathbf{Y} \mathbf{R}\_{i}(\mathbf{x}, \mathbf{k}) \right\| \tag{19}$$

time components of partial output quantities { i i ( )}i Y x ,k are obtained, in block SS1 reference quantities { i i ( )}i W x,∞ are computed. Then on the output of the algebraic block is { ii ii ii () ( )()}i E x ,k W x , Y x ,k = ∞− . These sequences { i i ( )}i E x ,k enter into the TS on{ i ( )}i R z and give { i ( )}i U k , which enter into HLDS with { i i ( )}i Y x ,k on the SS2 output - among { i ( )}i U k and { i i ( )}i Y x ,k there are relations { i i ( )}i SH x ,z . - This analysis of control synthesis process shows that synthesis in time domain is realized on the level of one parameter control loops { ii i ( ) ()}i SH x ,z ,R z , Fig. 9.

Fig. 11. Distributed parameter closed-loop system of control with reduced space components of output quantity

Control of Distributed Parameter Systems - Engineering Methods

Fig. 13. Lumped parameter SISO control loops – i-th control loop

SH x ,z - i-th time component of HLDS dynamics i i ( )

Y x ,k /U k ii i ( ) () - i-th controlled/control quantity

W x ,k /E x ,k ii ii ()() - i-th desired quantity/control error

equation (9), is accomplished as given by relation (16).

*Space Tasks of Control Synthesis* – as approximation tasks.

R z - i-th lumped parameter controller i ( )

synthesis is realized as:

in timestep k is obtained

the block SS2 -{ i i ( )}i Y x ,k

where E x,k ( )

approximation set { i ( )}i,k YR x,k

and Software Support in the MATLAB & Simulink Programming Environment 41

If for k → ∞ { i i () ( )}i,k YR x,k HR x, → ∞ *H* ,{ ii ii ()()} { ii ii () ()} i,k i Y x ,k Y x ,k Y x , W x , = → ∞= ∞

along with { ii ii () ( ) }i,k E x ,k E x , 0 → ∞= - this actually means that the control task defined in

Finally we may state as a summary, that in closed-loop control with RHLDS the control

At the same time the solution of the approximation problem in block SS1 on the

i 1 min W x, W x ,k YR x,k =

( ) ( ) () n ii i

( ) ( ) ()() n ii i

is the unique element at the best approximation of W x, ( ) ∞ on the set of

in the timestep k distributed output quantity Y x,k is given as ( )

W x, W x ,k YR x,k E x,k

approximate functions { i ( )}i YR x,k . Similary by the solution of approximation problem in

() ( ) () n ii i

i 1 Y x,k Y x ,k YR x,k

=

∞ − (20)

= (22)

∞ = + (21)

*Time Tasks of Control Synthesis* – on the level of lumped parameter control loops

i 1

=

HLDS - controlled system with zero-order hold units RHLDS - model for reduced space components: { i ( )}i,k YR x,k CS - control synthesis TS/SS1,SS2 - time/space parts of control synthesis K - time/space sampling Y x,t ( ) - distributed output quantity W x, ,V x,t ( )( ) ∞ - reference and disturbance quantities { i i ( )} { i i ( )} i,k i,k Y x ,k Y x ,k = – time components of output quantity { i ii ( )}i W W x, = ∞ /{ i i ( )}i,k E x ,k - reference quantities/control errors U k( ) - vector of lumped control quantities

Fig. 12. The block of time synthesis

TS - time part of control synthesis

{ i ( )}i R z - lumped parameter controllers

() ( ) { i i }i,k E k E x ,k = - vector of lumped control errors

() () { i }i,k Uk U k = - vector of lumped control quantities

For k → ∞ { i i () ( )}i,k YR x,k HR x, → ∞ *H* , { i i ( )} { ii ii () ()} i,k i Y x ,k Y x , W x , → ∞= ∞ along with { ii ii () ( ) }i,k E x ,k E x , 0 → ∞= . Thus the control task, defined by equation (9) is accomplished as given by relation (16). In case of the uncertainty of the control process relations similar to (17-18) are also valid.

Let´s consider now the approximation of W x, ( ) ∞ in the block SS1 in timestep k, on the set of { i ( )}i,k YR x,k . Then in the control process sequences of quantities { i i ( )}i,k W x ,k are obtained, as desired quantities of SISO control loops { ii i ( ) ()}i SH x ,z ,R z which are closed throughout the blocks TS, HLDS and SS2, see Fig. 13.

Fig. 13. Lumped parameter SISO control loops – i-th control loop

SH x ,z - i-th time component of HLDS dynamics i i ( )

R z - i-th lumped parameter controller i ( )

40 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

For k → ∞ { i i () ( )}i,k YR x,k HR x, → ∞ *H* , { i i ( )} { ii ii () ()} i,k i Y x ,k Y x , W x , → ∞= ∞ along with

{ ii ii () ( ) }i,k E x ,k E x , 0 → ∞= . Thus the control task, defined by equation (9) is accomplished as given by relation (16). In case of the uncertainty of the control process relations similar to

Let´s consider now the approximation of W x, ( ) ∞ in the block SS1 in timestep k, on the set of { i ( )}i,k YR x,k . Then in the control process sequences of quantities { i i ( )}i,k

obtained, as desired quantities of SISO control loops { ii i ( ) ()}i SH x ,z ,R z which are closed

W x ,k are

HLDS - controlled system with zero-order hold units

TS/SS1,SS2 - time/space parts of control synthesis

W x, ,V x,t ( )( ) ∞ - reference and disturbance quantities

{ i i ( )} { i i ( )} i,k i,k Y x ,k Y x ,k = – time components of output quantity

CS - control synthesis

K - time/space sampling

Y x,t ( ) - distributed output quantity

U k( ) - vector of lumped control quantities

{ i ii ( )}i W W x, = ∞ /{ i i ( )}i,k E x ,k

Fig. 12. The block of time synthesis TS - time part of control synthesis

() ( ) { i i }i,k E k E x ,k =

(17-18) are also valid.

{ i ( )}i R z - lumped parameter controllers

() () { i }i,k Uk U k = - vector of lumped control quantities

throughout the blocks TS, HLDS and SS2, see Fig. 13.

RHLDS - model for reduced space components: { i ( )}i,k YR x,k

Y x ,k /U k ii i ( ) () - i-th controlled/control quantity

W x ,k /E x ,k ii ii ()() - i-th desired quantity/control error

$$\begin{aligned} \text{1\text{\textquotedblleft for}} \quad \mathbf{k} \rightarrow \approx \left[ \mathbf{Y} \mathbf{R}\_{\mathbf{i}} (\mathbf{x}, \mathbf{k}) \rightarrow \mathcal{M} \mathbf{H} \mathbf{R}\_{\mathbf{i}} (\mathbf{x}\_{\circ} \boldsymbol{\approx}) \right]\_{\mathbf{i}, \mathbf{k} \; \prime} \left[ \mathbf{Y}\_{\mathbf{i}} (\mathbf{x}\_{\circ}, \mathbf{k}) = \mathbf{\bar{Y}}\_{\mathbf{i}} (\mathbf{x}\_{\circ}, \mathbf{k}) \right]\_{\mathbf{i}, \mathbf{k}} \rightarrow \left[ \mathbf{\bar{Y}}\_{\mathbf{i}} (\mathbf{x}\_{\circ}, \boldsymbol{\approx}) = \mathbf{\bar{W}}\_{\mathbf{i}} (\mathbf{x}\_{\circ}, \boldsymbol{\approx}) \right]\_{\mathbf{i}, \mathbf{k}} \end{aligned}$$

along with { ii ii () ( ) }i,k E x ,k E x , 0 → ∞= - this actually means that the control task defined in equation (9), is accomplished as given by relation (16).

Finally we may state as a summary, that in closed-loop control with RHLDS the control synthesis is realized as:

*Time Tasks of Control Synthesis* – on the level of lumped parameter control loops *Space Tasks of Control Synthesis* – as approximation tasks.

At the same time the solution of the approximation problem in block SS1 on the approximation set { i ( )}i,k YR x,k

$$\min \left\| \mathbf{W} \left( \mathbf{x}, \boldsymbol{\omega} \right) - \sum\_{i=1}^{n} \mathbf{W}\_{i} \left( \mathbf{x}\_{i}, \mathbf{k} \right) \mathbf{Y} \mathbf{R}\_{i} \left( \mathbf{x}, \mathbf{k} \right) \right\|\tag{20}$$

in timestep k is obtained

$$\mathbf{W}(\mathbf{x},\boldsymbol{\omega}\circ\mathbf{s})=\sum\_{i=1}^{n}\bar{\mathbf{W}}\_{i}(\mathbf{x}\_{i},\mathbf{k})\mathbf{Y}\mathbf{R}\_{i}(\mathbf{x},\mathbf{k})+\bar{\mathbf{E}}(\mathbf{x},\mathbf{k})\tag{21}$$

where E x,k ( ) is the unique element at the best approximation of W x, ( ) ∞ on the set of approximate functions { i ( )}i YR x,k . Similary by the solution of approximation problem in the block SS2 -{ i i ( )}i Y x ,k in the timestep k distributed output quantity Y x,k is given as ( )

$$\mathbf{Y}(\mathbf{x}, \mathbf{k}) = \sum\_{i=1}^{n} \bar{\mathbf{Y}}\_{i}(\mathbf{x}\_{i}, \mathbf{k}) \mathbf{Y} \mathbf{R}\_{i}(\mathbf{x}, \mathbf{k}) \tag{22}$$

#### **4.3 Closed-loop control**

Let us now consider a distributed parameter feedback loop as featured in Fig. 14. with initial conditions as above, where in timestep k an approximation problem is solved

$$\min \left\| \mathbb{E} \left( \mathbf{x}, \mathbf{k} \right) - \sum\_{i=1}^{n} \mathbb{E}\_{i} \left( \mathbf{x}\_{i}, \mathbf{k} \right) \mathbf{Y} \mathbf{R}\_{i} \left( \mathbf{x}, \mathbf{k} \right) \right\|\tag{23}$$

Control of Distributed Parameter Systems - Engineering Methods

Y x,t /W x, () ( ) ∞ – distributed controlled/reference quantity

{ i ( )}i YR x,k /{ i ( )}i *H*HR x,∞ - reduced characteristics in space domain

( ) ()() ( ) ( ) ()() n n

W x ,k YR x,k E x,k Y x ,k E x ,k YR x,k E x,k

{ ii ii ii ()()()}i,k W x ,k Y x ,k E x ,k = + is obtained. Now in vector form this means

ii i ii ii i

( ) () ( ) ( ) () n n ii i ii ii i

= +

W x ,k YR x,k Y x ,k E x ,k YR x,k

+= + +

CS – control synthesis

E k( )

{ i ( )}i

K – time/space sampling

TS /SS – time/space control synthesis

E x,k ( ) – distributed control error

– vector of control errors U k( ) – vector of control quantities

i 1 i 1

= =

i 1 i 1

Then sequences { i i ( )}i,k E x ,k

U k and { i i ( )}i Y x ,k

= =

V x,t ( ) – distributed disturbance quantity

and Software Support in the MATLAB & Simulink Programming Environment 43

Comparison of relation (21) and result of the approximation problem (24) gives

and then

, finally

( ) ( ) ( ) ( ) ( ) () ( ) ( ) ( ) { ii ii ii }i,k W k Y k E k E k W k Y k E x ,k W x ,k Y x ,k = + = −= = − (25)

synthesis process shows that synthesis in time domain is realized on the level of one parameter control loops { ii i ( ) ()}i SH x ,z ,R z , Fig. 13. - closed throughout the structure of distributed parameter control loop, Fig. 14. If for k → ∞ { i i () ( )}i YR x,k HR x, → ∞ *H* and

( ) ( ) () ( )( ) ( ) ()

ii i ii i ii i

W x , Y x , HR x, E x, E x , HR x, E x,

= ∞ ∞ ∞+ ∞− ∞ ∞ = ∞

W x , HR x, E x, Y x , HR x, E x , HR x,

*H HH*

= ∞ ∞+ ∞− ∞ ∞− ∞ ∞ =

( ) ()( ) ( ) () ( ) ()

ii i ii i

*H H*

∞− ∞ ∞ = ∞− ∞− ∞ ∞ =

( )( ) ()( ) ( ) () ( )

is valid, so the above given control task (9) is accomplished - whereas in the steady-state { ii ii ii ( )( ) ( ) }i W x, Y x, E x, 0 ∞− ∞= ∞= . By concluding the above presented discussion, the

{ i i ( )} { ii ii () ()} i,k i Y x ,k Y x , W x , → ∞= ∞ along with { ii ii () ( ) }i,k E x ,k E x , 0 → ∞= then

min E x, E x , HR x, W x, Y x, E x , HR x,

n n

i 1 i 1

= =

ii ii i ii i

n n n

i 1 i 1 i 1

= = =

n n

i 1 i 1

= =

control synthesis in closed-loop control is realized as:

enter into the block TS on { i ( )}i R z and give { i ( )}i U k , among

there are relations { i i ( )}i SH x ,z . Finally this analysis of control

(26)

{ i i ( )}i SH x ,z - transfer functions - dynamic characteristics of HLDS in time domain

and as a result in timestep k a vector () ( ) { i i }i E k E x ,k = is obtained. By relations (20-22) the further equations are valid

$$\begin{aligned} &\min \left\| \mathbb{E} \left( \mathbf{x}, \mathbf{k} \right) - \sum\_{i=1}^{n} \mathbb{E}\_{i} \left( \mathbf{x}\_{i}, \mathbf{k} \right) \mathbf{Y} \mathbf{R}\_{i} \left( \mathbf{x}, \mathbf{k} \right) \right\| = \\ &\mathbf{x} = \min \left\| \mathbf{W} \left( \mathbf{x}, \boldsymbol{\omega} \right) - \mathbf{Y} \left( \mathbf{x}, \mathbf{k} \right) - \sum\_{i=1}^{n} \mathbb{E}\_{i} \left( \mathbf{x}\_{i}, \mathbf{k} \right) \mathbf{Y} \mathbf{R}\_{i} \left( \mathbf{x}, \mathbf{k} \right) \right\| = \\ &\mathbf{=} \min \left\| \mathbf{W} \left( \mathbf{x}, \boldsymbol{\omega} \right) - \sum\_{i=1}^{n} \left[ \bar{\mathbf{Y}}\_{i} \left( \mathbf{x}\_{i}, \mathbf{k} \right) + \bar{\mathbf{E}\_{i} \left( \mathbf{x}\_{i}, \mathbf{k} \right)} \right] \mathbf{Y} \mathbf{R}\_{i} \left( \mathbf{x}, \mathbf{k} \right) \right\| = \\ &\mathbf{=} \left\| \mathbf{W} \left( \mathbf{x}, \boldsymbol{\omega} \right) - \sum\_{i=1}^{n} \left[ \bar{\mathbf{Y}}\_{i} \left( \mathbf{x}\_{i}, \mathbf{k} \right) + \bar{\mathbf{E}\_{i} \left( \mathbf{x}\_{i}, \mathbf{k} \right)} \right] \mathbf{Y} \mathbf{R}\_{i} \left( \mathbf{x}, \mathbf{k} \right) \right\| \end{aligned} \tag{24}$$

The problem solution ( ) ( ) ( ) ()() n ii ii i i 1 W x, Y x ,k E x ,k YR x,k E x,k = ∞ = + + is obtained by approximation.

Fig. 14. Distributed parameter closed-loop system of control

HLDS - LDS with zero-order hold units

Let us now consider a distributed parameter feedback loop as featured in Fig. 14. with initial

i 1 min E x,k E x ,k YR x,k =

() ( ) ()

n

i 1

=

= ∞− +

The problem solution ( ) ( ) ( ) ()() n

i 1

=

Fig. 14. Distributed parameter closed-loop system of control

HLDS - LDS with zero-order hold units

ii i

− =

min W x, Y x,k E x ,k YR x,k

= ∞− − =

n

min E x,k E x ,k YR x,k

i 1

=

n

i 1

=

() ( ) () n ii i

( )() ( ) ()

n

i 1

=

( ) ( ) ( ) ()

W x, Y x ,k E x ,k YR x,k

min W x, Y x ,k E x ,k YR x,k

= ∞− + =

( ) ( ) ( ) ()

ii ii i

ii ii i

W x, Y x ,k E x ,k YR x,k E x,k

∞ = + +

ii ii i

ii i

− (23)

(24)

is obtained. By relations (20-22) the

is obtained by

conditions as above, where in timestep k an approximation problem is solved

and as a result in timestep k a vector () ( ) { i i }i E k E x ,k =

**4.3 Closed-loop control**

further equations are valid

approximation.

TS /SS – time/space control synthesis

K – time/space sampling

Y x,t /W x, () ( ) ∞ – distributed controlled/reference quantity

{ i i ( )}i SH x ,z - transfer functions - dynamic characteristics of HLDS in time domain

{ i ( )}i YR x,k /{ i ( )}i *H*HR x,∞ - reduced characteristics in space domain

E x,k ( ) – distributed control error

V x,t ( ) – distributed disturbance quantity

E k( ) – vector of control errors

U k( ) – vector of control quantities

Comparison of relation (21) and result of the approximation problem (24) gives ( ) ()() ( ) ( ) ()() n n ii i ii ii i i 1 i 1 W x ,k YR x,k E x,k Y x ,k E x ,k YR x,k E x,k = = += + + and then

$$\sum\_{i=1}^{n} \bar{\mathbf{W}}\_{i}(\mathbf{x}\_{i}, \mathbf{k}) \mathbf{Y} \mathbf{R}\_{i}(\mathbf{x}, \mathbf{k}) = \sum\_{i=1}^{n} \left[ \bar{\mathbf{Y}}\_{i}(\mathbf{x}\_{i}, \mathbf{k}) + \bar{\mathbf{E}}\_{i}(\mathbf{x}\_{i}, \mathbf{k}) \right] \mathbf{Y} \mathbf{R}\_{i}(\mathbf{x}, \mathbf{k}) \, , \tag{11.8}$$

$$\left\{ \bar{\mathbf{W}}\_{i}(\mathbf{x}\_{i}, \mathbf{k}) = \bar{\mathbf{Y}}\_{i}(\mathbf{x}\_{i}, \mathbf{k}) + \bar{\mathbf{E}}\_{i}(\mathbf{x}\_{i}, \mathbf{k}) \right\}\_{i, \mathbf{k}} \text{ is obtained. Now in vector form this means$$

$$\bar{\tilde{\mathbf{W}}}(\mathbf{k}) = \bar{\tilde{\mathbf{Y}}}(\mathbf{k}) + \bar{\tilde{\mathbf{E}}}(\mathbf{k}) \Rightarrow \bar{\tilde{\mathbf{E}}}(\mathbf{k}) = \bar{\tilde{\mathbf{W}}}(\mathbf{k}) - \bar{\tilde{\mathbf{Y}}}(\mathbf{k}) = \left[ \bar{\mathbf{E}}\_{\mathrm{i}}(\mathbf{x}\_{i}, \mathbf{k}) = \bar{\mathbf{W}}\_{\mathrm{i}}(\mathbf{x}\_{i}, \mathbf{k}) - \bar{\mathbf{Y}}\_{\mathrm{i}}(\mathbf{x}\_{i}, \mathbf{k}) \right]\_{\mathrm{i}, \mathbf{k}} \tag{25}$$

Then sequences { i i ( )}i,k E x ,k enter into the block TS on { i ( )}i R z and give { i ( )}i U k , among { i ( )}i U k and { i i ( )}i Y x ,k there are relations { i i ( )}i SH x ,z . Finally this analysis of control synthesis process shows that synthesis in time domain is realized on the level of one parameter control loops { ii i ( ) ()}i SH x ,z ,R z , Fig. 13. - closed throughout the structure of distributed parameter control loop, Fig. 14. If for k → ∞ { i i () ( )}i YR x,k HR x, → ∞ *H* and { i i ( )} { ii ii () ()} i,k i Y x ,k Y x , W x , → ∞= ∞ along with { ii ii () ( ) }i,k E x ,k E x , 0 → ∞= then

( ) ( ) () ( )( ) ( ) () ( ) ()( ) ( ) () ( ) () ( )( ) ()( ) ( ) () ( ) n n ii i ii i i 1 i 1 n n n ii i ii i ii i i 1 i 1 i 1 n n ii ii i ii i i 1 i 1 min E x, E x , HR x, W x, Y x, E x , HR x, W x , HR x, E x, Y x , HR x, E x , HR x, W x , Y x , HR x, E x, E x , HR x, E x, = = = = = = = ∞− ∞ ∞ = ∞− ∞− ∞ ∞ = = ∞ ∞+ ∞− ∞ ∞− ∞ ∞ = = ∞ ∞ ∞+ ∞− ∞ ∞ = ∞ *H H H HH* - *H H* (26)

is valid, so the above given control task (9) is accomplished - whereas in the steady-state { ii ii ii ( )( ) ( ) }i W x, Y x, E x, 0 ∞− ∞= ∞= . By concluding the above presented discussion, the control synthesis in closed-loop control is realized as:

Control of Distributed Parameter Systems - Engineering Methods

operation.

Corporation

and Software Support in the MATLAB & Simulink Programming Environment 45

The block **Tutorial** presents methodological framework for formulation and solution of control problems for DPS. The block **Show** contains motivation examples such as: *Control of temperature field of 3D metal body* (the controlled system was modelled in the virtual software environment COMSOL Multiphysics); *Control of 3D beam of "smart" structure* (the controlled system was modelled in the virtual software environment ANSYS)*; Adaptive control of glass furnace* (the controlled system was modelled by Partial Differential Equations Toolbox of the MATLAB ), and *Groundwater remediation control* (the controlled system was modelled in the virtual software environment MODFLOW)**.** The block **Demos** contains examples oriented at the methodology of modelling and control synthesis. The **DPS Wizard** gives an automatized guide for arrangement and setting distributed parameter control loops in step-by-step

Fig. 15. Distributed Parameter Systems Blockset on the web portal of The MathWorks

*Time Tasks of Control Synthesis* – on the level of lumped parameter control loops *Space Tasks of Control Synthesis* – as approximation tasks.

When mathematical models cannot provide an ideal description of controlled DPS dynamics and disturbances are present with an overall effect on the output in steady-state, expressed by EY x, ( ) ∞ then the realtions similar to (17-18) are also valid here.

In practice mostly only reduced distributed parameter transient responses in steady-state { i ( )}i *H*HR x,∞ are considered for the solution of the approximation tasks in the block SS of

the scheme in Fig. 14. along with robustification of controllers { i ( )}i Rz .

For simplicity problems of DPS control have been formulated here for the distributed desired quantity W x, ( ) ∞ . In case of W x,k is assumed, the control synthesis is realized ( ) similarly:

The solution of the presented problems of control synthesis an assumption is used, that in the framework of the chosen control systems the prescribed control quality can be reached both in the in space and time domain. However in the design of actual control systems for the given distributed parameter systems, usually the

is required and necessary.

#### **5. Distributed Parameter Systems Blockset for MATLAB & Simulink**

As a software support for DPS modelling, control and design of problems in MATLAB & Simulink the programming environment **Distributed Parameter Systems Blockset for MATLAB & Simulink** (**DPS Blockset**) - a Third-Party Product of The MathWorks www.mathworks.com/products/connections/ – Fig. 15., has been developed within the program CONNECTIONS of The MathWorks Corporation by the Institute of Automation, Measurement and Applied Informatics of Mechanical Engineering Faculty, Slovak University of Technology in Bratislava (IAMAI-MEF-STU) (Hulkó et al., 2003- 2010). Fig. 16. shows The library of **DPS Blockset**. The **HLDS** and **RHLDS** blocks model controlled DPS dynamics described by numerical structures as LDS with zero-order hold units - H. **DPS Control Synthesis** provides feedback to distributed parameter controlled systems in control loops with blocks for discrete-time **PID, Algebraic, State-Space and Robust Synthesis**. The block **DPS Input** generates distributed quantities, which can be used as distributed reference quantities or distributed disturbances, etc. **DPS Display** presents distributed quantities with many options including export to AVI files. The block **DPS Space Synthesis** performs space synthesis as an approximation problem.

As a demonstration, some results of the discrete-time PID control of complex-shape metal body heating by the **DPS Blockset** are shown in Fig. 17.-19., where the heating process was modelled by finite element method in the COMSOL Multiphysics virtual software environment - www.comsol.com.

When mathematical models cannot provide an ideal description of controlled DPS dynamics and disturbances are present with an overall effect on the output in steady-state, expressed

In practice mostly only reduced distributed parameter transient responses in steady-state { i ( )}i *H*HR x,∞ are considered for the solution of the approximation tasks in the block SS of

For simplicity problems of DPS control have been formulated here for the distributed desired quantity W x, ( ) ∞ . In case of W x,k is assumed, the control synthesis is realized ( )

The solution of the presented problems of control synthesis an assumption is used, that in the framework of the chosen control systems the prescribed control quality can be reached both in the in space and time domain. However in the design of actual control systems for

• optimization of dynamical characteristics of lumped/distributed parameter actuators

As a software support for DPS modelling, control and design of problems in MATLAB & Simulink the programming environment **Distributed Parameter Systems Blockset for MATLAB & Simulink** (**DPS Blockset**) - a Third-Party Product of The MathWorks www.mathworks.com/products/connections/ – Fig. 15., has been developed within the program CONNECTIONS of The MathWorks Corporation by the Institute of Automation, Measurement and Applied Informatics of Mechanical Engineering Faculty, Slovak University of Technology in Bratislava (IAMAI-MEF-STU) (Hulkó et al., 2003- 2010). Fig. 16. shows The library of **DPS Blockset**. The **HLDS** and **RHLDS** blocks model controlled DPS dynamics described by numerical structures as LDS with zero-order hold units - H. **DPS Control Synthesis** provides feedback to distributed parameter controlled systems in control loops with blocks for discrete-time **PID, Algebraic, State-Space and Robust Synthesis**. The block **DPS Input** generates distributed quantities, which can be used as distributed reference quantities or distributed disturbances, etc. **DPS Display** presents distributed quantities with many options including export to AVI files. The block **DPS Space Synthesis** performs space synthesis as an approximation

As a demonstration, some results of the discrete-time PID control of complex-shape metal body heating by the **DPS Blockset** are shown in Fig. 17.-19., where the heating process was modelled by finite element method in the COMSOL Multiphysics virtual software

• optimization of dynamical characteristics of lumped parameter control loops

**5. Distributed Parameter Systems Blockset for MATLAB & Simulink**

Rz .

*Time Tasks of Control Synthesis* – on the level of lumped parameter control loops

*Space Tasks of Control Synthesis* – as approximation tasks.

similarly:

problem.

environment - www.comsol.com.

is required and necessary.

by EY x, ( ) ∞ then the realtions similar to (17-18) are also valid here.

the scheme in Fig. 14. along with robustification of controllers { i ( )}i

throughout structures of the distributed parameter control loop.

the given distributed parameter systems, usually the • optimization of the number and layout of actuators The block **Tutorial** presents methodological framework for formulation and solution of control problems for DPS. The block **Show** contains motivation examples such as: *Control of temperature field of 3D metal body* (the controlled system was modelled in the virtual software environment COMSOL Multiphysics); *Control of 3D beam of "smart" structure* (the controlled system was modelled in the virtual software environment ANSYS)*; Adaptive control of glass furnace* (the controlled system was modelled by Partial Differential Equations Toolbox of the MATLAB ), and *Groundwater remediation control* (the controlled system was modelled in the virtual software environment MODFLOW)**.** The block **Demos** contains examples oriented at the methodology of modelling and control synthesis. The **DPS Wizard** gives an automatized guide for arrangement and setting distributed parameter control loops in step-by-step operation.

Fig. 15. Distributed Parameter Systems Blockset on the web portal of The MathWorks Corporation

Control of Distributed Parameter Systems - Engineering Methods

metal body in DPS Blockset environment

numerical net

environment

and Software Support in the MATLAB & Simulink Programming Environment 47

Fig. 17. Distributed parameter control loop for discrete-time PID control of heating of

Fig. 18. Distributed reference and controlled quantities of metal body heating over the

Fig. 19. Quadratic norm of distributed control error and discrete lumped actuating quantities at discrete-time PID control of heating of metal body in DPS Blockset

Fig. 16. The library of Distributed Parameter Systems Blockset for MATLAB & Simulink – Third-Party Product of The MathWorks

Fig. 16. The library of Distributed Parameter Systems Blockset for MATLAB & Simulink –

Third-Party Product of The MathWorks

Fig. 17. Distributed parameter control loop for discrete-time PID control of heating of metal body in DPS Blockset environment

Fig. 18. Distributed reference and controlled quantities of metal body heating over the numerical net

Fig. 19. Quadratic norm of distributed control error and discrete lumped actuating quantities at discrete-time PID control of heating of metal body in DPS Blockset environment

Control of Distributed Parameter Systems - Engineering Methods

means of the Optimization Toolbox.

DPS via the Internet.

**9. References**

**8. Acknowledgment**

*controlled distributed parameter systems".*

*of 8-th World Congress of IFAC*, Kyoto, 1981

(in Russian)

Munich, 1987

*Systems*, Perpignan, 1989

and Software Support in the MATLAB & Simulink Programming Environment 49

for the analogous decomposition of control synthesis into the space and time subtasks. In space domain approximation problems are solved, while in the time domain control synthesis is realized by lumped parameter SISO control loops (Hulkó et al., 1981-2010). Based on these decomposition a software product named **Distributed Parameter Systems Blockset for MATLAB & Simulink** - a Third-Party software product of The MathWorks www.mathworks.com/products/connections/ has been developed within the program CONNECTIONS of The MathWorks Corporation, (Hulkó et al., 2003-2010), where time domain toolboxes and blocksets of software environment MATLAB & Simulink as Control Systems Toolbox, Simulink Control Design, System Identification Toolbox,... are made use of. In the space domain approximation problems are solved as optimization problems by

For the further support of research in this area a web portal named **Distributed Parameter Systems Control** - www.dpscontrol.sk was realized (Hulkó et al., 2003-2010), see Fig. 20. for an illustration. On the above mentioned web portal, the online version of the monograph titled **Modeling, Control and Design of Distributed Parameter Systems with Demonstrations in MATLAB** - www.mathworks.com/support/books/ (Hulkó et al., 1998), is presented along with application examples from different disciplines such as: control of technological and production processes, control and design of mechatronic structures, groundwater remediation control, etc. This web portal also offers for those interested the download of the demo version of the **Distributed Parameter Systems Blockset for MATLAB & Simulink** with **Tutorial** , **Show** , **Demos** and **DPS Wizard**. This portal also offers the **Interactive Control** service for interactive solution of model control problems of

This work was supported by the Slovak Scientific Grant Agency VEGA under the contract No. 1/0138/11 for project *"Control of dynamical systems given by numerical structures as distributed parameter systems"* and the Slovak Research and Development Agency under the contract No. APVV-0160-07 for project "*Advanced Methods for Modeling, Control and Design of Mechatronical Systems as Lumped-input and Distributed-output Systems"* also the project No*.*  APVV-0131-10 *"High-tech solutions for technological processes and mechatronic components as*

Butkovskij, A. G. (1965). *Optimal control of distributed parameter systems*. Nauka, Moscow

Hulkó, G. et al. (1981). On Adaptive Control of Distributed Parameter Systems, *Proceedings*

Hulkó, G. et al. (1987). Control of Distributed Parameter Systems by means of Multi-Input

Hulkó, G. (1989). Identification of Lumped Input and Distributed Output Systems,

and Multi-Distributed-Output Systems*, Proceedings of 10-th World Congress of IFAC*,

*Proceedings of 5-th IFAC / IMACS / IFIP Symposium on Control of Distributed Parameter*

#### **6. Interactive control via the Internet**

For the interactive formulation and solution of DPS demonstration control problems via the Internet, an **Interactive Control** service has been started on the web portal **Distributed Parameter Systems Control -** www.dpscontrol.sk of the IAMAI-MEF-STU (Hulkó, 2003-2010) – see Fig. 20. for a screenshot of the site. In the framework of the problem formulation, first the computational geometry and mesh are chosen in the complex 3D shape definition domain, then DT distributed transient responses are computed in virtual software environments for numerical dynamical analysis of machines and processes. Finally, the distributed reference quantity is specified in points of the computational mesh - Fig. 18. Representing the solution to those interested animated results of actuating quantities, quadratic norm of control error, distributed reference and controlled quantity are sent in the form of **DPS Blockset** outputs – see Fig. 17-19. for illustration.

Fig. 20. Web portal Distributed Parameter Systems Control with monograph Modeling, Control and Design of Distributed Parameter Systems with Demonstrations in MATLAB and service Interactive Control

#### **7. Conclusion**

The aim of this chapter is to present a philosophy of the engineering approach for the control of DPS – given by numerical structures, which opens a wide space for novel applications of the toolboxes and blocksets of the MATLAB & Simulink software environment. This approach is based on the general decomposition into time and space components of controlled DPS dynamics represented by numerically computed distributed parameter transient and impulse characteristics, given on complex shape definition domains in 3D. Starting out from this dynamics decomposition a methodical framework is presented for the analogous decomposition of control synthesis into the space and time subtasks. In space domain approximation problems are solved, while in the time domain control synthesis is realized by lumped parameter SISO control loops (Hulkó et al., 1981-2010). Based on these decomposition a software product named **Distributed Parameter Systems Blockset for MATLAB & Simulink** - a Third-Party software product of The MathWorks www.mathworks.com/products/connections/ has been developed within the program CONNECTIONS of The MathWorks Corporation, (Hulkó et al., 2003-2010), where time domain toolboxes and blocksets of software environment MATLAB & Simulink as Control Systems Toolbox, Simulink Control Design, System Identification Toolbox,... are made use of. In the space domain approximation problems are solved as optimization problems by means of the Optimization Toolbox.

For the further support of research in this area a web portal named **Distributed Parameter Systems Control** - www.dpscontrol.sk was realized (Hulkó et al., 2003-2010), see Fig. 20. for an illustration. On the above mentioned web portal, the online version of the monograph titled **Modeling, Control and Design of Distributed Parameter Systems with Demonstrations in MATLAB** - www.mathworks.com/support/books/ (Hulkó et al., 1998), is presented along with application examples from different disciplines such as: control of technological and production processes, control and design of mechatronic structures, groundwater remediation control, etc. This web portal also offers for those interested the download of the demo version of the **Distributed Parameter Systems Blockset for MATLAB & Simulink** with **Tutorial** , **Show** , **Demos** and **DPS Wizard**. This portal also offers the **Interactive Control** service for interactive solution of model control problems of DPS via the Internet.

### **8. Acknowledgment**

48 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

For the interactive formulation and solution of DPS demonstration control problems via the Internet, an **Interactive Control** service has been started on the web portal **Distributed Parameter Systems Control -** www.dpscontrol.sk of the IAMAI-MEF-STU (Hulkó, 2003-2010) – see Fig. 20. for a screenshot of the site. In the framework of the problem formulation, first the computational geometry and mesh are chosen in the complex 3D shape definition domain, then DT distributed transient responses are computed in virtual software environments for numerical dynamical analysis of machines and processes. Finally, the distributed reference quantity is specified in points of the computational mesh - Fig. 18. Representing the solution to those interested animated results of actuating quantities, quadratic norm of control error, distributed reference and controlled quantity are

sent in the form of **DPS Blockset** outputs – see Fig. 17-19. for illustration.

Fig. 20. Web portal Distributed Parameter Systems Control with monograph Modeling, Control and Design of Distributed Parameter Systems with Demonstrations in MATLAB

The aim of this chapter is to present a philosophy of the engineering approach for the control of DPS – given by numerical structures, which opens a wide space for novel applications of the toolboxes and blocksets of the MATLAB & Simulink software environment. This approach is based on the general decomposition into time and space components of controlled DPS dynamics represented by numerically computed distributed parameter transient and impulse characteristics, given on complex shape definition domains in 3D. Starting out from this dynamics decomposition a methodical framework is presented

**6. Interactive control via the Internet**

and service Interactive Control

**7. Conclusion**

This work was supported by the Slovak Scientific Grant Agency VEGA under the contract No. 1/0138/11 for project *"Control of dynamical systems given by numerical structures as distributed parameter systems"* and the Slovak Research and Development Agency under the contract No. APVV-0160-07 for project "*Advanced Methods for Modeling, Control and Design of Mechatronical Systems as Lumped-input and Distributed-output Systems"* also the project No*.*  APVV-0131-10 *"High-tech solutions for technological processes and mechatronic components as controlled distributed parameter systems".*

#### **9. References**

**3**

Lubomír Brančík

*Czech Republic*

α*tT*),

}, with **c** = (*c*1,...,*cn*) as an

*Brno University of Technology*

**Numerical Inverse Laplace Transforms for**

Numerical inverse Laplace transform (NILT) methods are widely used in various scientific areas, especially for a solution of respective differential equations. In field of an electrical engineering many various approaches have been considered so far, but mostly for a single variable (1D NILT), see at least (Brančík, 1999, 2007b; Cohen, 2007; Valsa & Brančík, 1998; Wu at al., 2001) from plenty of papers. Much less attention was paid to multidimensional variable (*n*D NILT) methods, see e.g. (Hwang at al., 1983; Singhal at al., 1975), useful rather for more complicated electromagnetic systems. The 2D NILT methods, see e.g. (Brančík, 2005, 2007a, 2007b; Hwang & Lu, 1999), can be applied for a transmission line analysis, or *n*D NILT methods, *n* ≥ 2, for a nonlinear circuits analysis, if relevant Laplace transforms are developed through a Volterra series expansion, see e.g. (Brančík, 2010a, 2010b, Karmakar, 1980; Schetzen, 2006), to highlight at least a few applications. This paper is focused on the class of NILT methods based on complex Fourier series approximation, their error analysis, their effective algorithms development in a Matlab language, and after all, on their selected applications in field of electrical engineering to show practical usefulness of the algorithms.

An *n*-dimensional Laplace transform of a real function *f*(*t*), with *t* = (*t*1,...,*tn*) as a row vector

0 0 1

abscissa of convergence, and the inequality taken componentwise, the original function is

*i n- fold i F f dt*

1 1

−∞ −∞ = ∏ *n*

*c n T*

*n T*

=

*s t st* − ∏ , (1)

*n*) being a minimal abscissa of convergence, and the *n*D

α

*i*

*i*

*F dt*

*t s st* . (2)

( ) = ( )exp( )

where *s* = (*s*1,...,*sn*) and *T* means a transposition. Under an assumption |*f*(*t*)| < *M*exp(

1 ( )= ( )exp( ) (2 )

*cj cj*

*n*

+ ∞ + ∞

*j c j*

∞ ∞

1

*n*

**2. Multidimensional numerical inverse Laplace transform**

of *n* real variables, is defined as (Hwang at al., 1983)

α = (α1,...,α

*f*

Laplace transform *F*(*s*) defined on a region {*s* ∈ *Cn*: Re[*s*] >

π

*j*

with *M* real positive and

given by an *n*-fold Bromwich integral

**1. Introduction**

**Electrical Engineering Simulation**

### **Numerical Inverse Laplace Transforms for Electrical Engineering Simulation**

Lubomír Brančík

*Brno University of Technology Czech Republic*

#### **1. Introduction**

50 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

Hulkó, G. et al. (1990). Computer Aided Design of Distributed Parameter Systems of

Hulkó G. (1991). Lumped Input and Distributed Ouptut Systems at the Control of

Hulkó, G. et al. (1998). *Modeling, Control and Design of Distributed Parameter Systems with*

Hulkó, G. et al. (2005). Web-based control design environment for distributed parameter

Hulkó, G. et al. (2003-2010). *Distributed Parameter Systems Control*. Web portal, Available

Hulkó, G. et al. (2003-2010). *Distributed Parameter Systems Blockset for MATLAB & Simulink*,

Hulkó, G. et al. (2009). Engineering Methods and Software Support for Modelling and

Hulkó, G. et. al (2010). Control of Technological Processes Modelled in COMSOL

Lasiecka, I., Triggiani, R. (2000). *Control Theory for Partial Differential Equations* (Encyclopedia of Mathematics and Its Applications 74), Cambridge U. Press, Cambridge, UK Lions, J. L. (1971). *Optimal control of systems governed by partial differential equations,* Springer-

Shadrin, A. (2005). *Approximation theory.* DAMTP University of Cambridge, Cambridge UK,

Wang, P. K. C. (1964). *Control of distributed parameter systems* (Advances in Control Systems:

Theory and Applications, 1.), Academic Press, New York

MathWorks, Bratislava-Natick, Available from: www.dpscontrol.sk Hulkó, G. et al. (2009). Engineering Methods and Software Support for Modeling and

*Proceedings of European Control Conference 2009*, Budapest, 2009

Distributed Parameter Systems. *Problems of Control and Information Theory,* Vol. 20,

*Demonstrations in MATLAB*, Publishing House STU, ISBN 80-227-1083-0, Bratislava

systems control education*, Proceedings of 16-th World Congress of IFAC*, Prague, 2005

www.mathworks.com/products/connections/ - Third-Party Product of The

Design of Discrete-time Control of Distributed Parameter Systems, *Mini-tutorial*,

Design of Discrete-time Control of Distributed Parameter Systems*. European Journal of Control,* Vol. 15, No. Iss. 3-4, *Fundamental Issues in Control,* (May-August 2009),

Multiphysics as Distributed Parameter Systems*, Proceedings of Asian COMSOL*

Control, *Proceedings of 11-th World Congress of IFAC,* Tallin, 1990

No. 2, pp. 113-128, Pergamon Press, Oxford

from: www.dpscontrol.sk

pp. 407-417, ISSN 0947-3580

*Conference,* Bangalore, 2010

Verlag, Berlin - Heidelberg - New York

Available from: www.damtp.cam.ac.uk

Numerical inverse Laplace transform (NILT) methods are widely used in various scientific areas, especially for a solution of respective differential equations. In field of an electrical engineering many various approaches have been considered so far, but mostly for a single variable (1D NILT), see at least (Brančík, 1999, 2007b; Cohen, 2007; Valsa & Brančík, 1998; Wu at al., 2001) from plenty of papers. Much less attention was paid to multidimensional variable (*n*D NILT) methods, see e.g. (Hwang at al., 1983; Singhal at al., 1975), useful rather for more complicated electromagnetic systems. The 2D NILT methods, see e.g. (Brančík, 2005, 2007a, 2007b; Hwang & Lu, 1999), can be applied for a transmission line analysis, or *n*D NILT methods, *n* ≥ 2, for a nonlinear circuits analysis, if relevant Laplace transforms are developed through a Volterra series expansion, see e.g. (Brančík, 2010a, 2010b, Karmakar, 1980; Schetzen, 2006), to highlight at least a few applications. This paper is focused on the class of NILT methods based on complex Fourier series approximation, their error analysis, their effective algorithms development in a Matlab language, and after all, on their selected applications in field of electrical engineering to show practical usefulness of the algorithms.

#### **2. Multidimensional numerical inverse Laplace transform**

An *n*-dimensional Laplace transform of a real function *f*(*t*), with *t* = (*t*1,...,*tn*) as a row vector of *n* real variables, is defined as (Hwang at al., 1983)

$$F(\mathbf{s}) = \bigcap\_{0}^{\infty} \dots \bigcap\_{n-\text{fold}}^{\infty} f(\mathbf{t}) \exp(-\mathbf{s}\mathbf{t}^T) \prod\_{i=1}^{n} dt\_i \,\,\tag{1}$$

where *s* = (*s*1,...,*sn*) and *T* means a transposition. Under an assumption |*f*(*t*)| < *M*exp(α*tT*), with *M* real positive and α = (α1,...,α*n*) being a minimal abscissa of convergence, and the *n*D Laplace transform *F*(*s*) defined on a region {*s* ∈ *Cn*: Re[*s*] > α}, with **c** = (*c*1,...,*cn*) as an abscissa of convergence, and the inequality taken componentwise, the original function is given by an *n*-fold Bromwich integral

$$f(t) = \frac{1}{\left(2\pi j\right)^{n}} \int\_{c\_{1} - j^{\infty}}^{c\_{1} + j^{\infty}} \cdots \int\_{c\_{n} - j^{\infty}}^{c\_{n} + j^{\infty}} F(\mathbf{s}) \exp(\mathbf{s}\mathbf{t}^{T}) \prod\_{i=1}^{n} dt\_{i} \,\,. \tag{2}$$

Numerical Inverse Laplace Transforms for Electrical Engineering Simulation 53

*n*

with *si* = *ci* + *jmi*Ω*i*. If we define *Fn* ≡ *F*(*s*1,...,*sn*-1,*sn*) and *F*0 ≡ *f*(*t*1,...,*tn*-1,*tn*), then *n* consequential

1 1 21 1

−− − −

{ } 1

1 1 ( )= ( ) *i i*

*c kT*

*i i i i i m e F F e*

τ

2 ( )= ( ) = ( )

*i i i ii ii i*

can be used *n* times, *i* = *n*,*n*-1...,1, to evaluate (8), with *pn* = (*s*1,...,*sn*-1,*sn*), *pn*-1 = (*s*1,...,*sn*-1,*tn*),..., *p*1 = (*s*1,...,*tn*-1,*tn*), and *p*0 = (*t*1,...,*tn*-1,*tn*), while *pj* = *sj* for *j* ≤ *i*, and *pj* = *tj* otherwise. A further technique is based on demand to find the solution on a whole region of discrete points. Then, taking into account *tik* = *kTi* in (9), with *Ti* as the sampling periods in the original

∞

=−∞

*i* = 1,...,*n*, a maximum reachable region is Οmax = [0,(*M*1-1)*T*1] ×...× [0,(*Mn*-1)*Tn*]. Thus, to meet

region of the calculation is chosen to be Οcal = [0,*t*1cal] ×...× [0,*tn*cal], with *ti*cal = (*Mi*/2–1)*Ti*,

As is shown in (Brančík, 2007a, 2010c), the discretized formula (10) can be evaluated by the FFT and IFFT algorithms, in conjunction with the quotient-difference (q-d) algorithm for

To explain it in more detail, let us consider an *r*-th cycle in gaining the original function via

*r r m*

+ . (11)

accelerating convergence of the residual infinite series, see following procedures.

− − *p p* . For its discretized version (10) we have

*c kT*

τ

2 ( , , , , )= ( , , , , ) *r r*

=−∞

*r rn r r n*

*e F s kT t F s c jm t e*

∞

*F F F e ds* π*j*

*nn n n n n n n nn*

*F Fs s t F Fs tt*

11 1 1

*F ft t t*

As is obvious we need to use a procedure able to make the inversion of Laplace transforms dependent on another *n-*1 parameters, complex in general. Let us denote arguments in (8) by

*ee e ft t t Fs s s e e e*

 τ

*c t c t c t jm t jm t jm t*

11 1

− −

= ( , , , ).

1

*n n*

−

*i*

*c j*

+ ∞

− ∞

2

*j mkT*

π τ

*i i*

*i* substituted. As follows from the error analysis (Brančík,

1) ×...× [0,

τ

2

*j mkT*

π τ

*r r*

*n*). For *k* = 0*,*1*,...,Mi*-1,

*i* = *MiTi*, *i* = 1,...,*n*. In practice, a

*p p* , (10)

τ

τ

π

 τ

*i*

*c j*

*i i*

*s t*

*pp p* (9)

= ( , , , ), = ( , , , ),

∞∞ ∞ Ω ΩΩ

2 22 1 11

(8)

*n nn*

, (7)

1 1 2 2

1 2

ττ

partial inversions are performed as

*pi* = (*p*1,...,*pn*-1,*pn*). Then the ILT of the type

domain, we can write an approximate formula

*i* = *n*,*n*-1...,1, and with Ω*i* = 2

(9), i.e. { } 1 1 1 ( )= ( ) *F F r r r rr* −

−

*i* = 1,...,*n*, to provide certain margins.

1 2 1 2 1 2

( , , , )= (,, ,) *n n*

*n n mm m n*

> { } { }

1

−

1

−

1 1

− −

π/τ

2011) a relative error is predictable on the region Οerr = [0,

the necessary condition Οmax ⊂ Οerr, we can set up fittingly

**2.2.2 FFT, IFFT, and quotient-difference algorithms utilization**

1 1 1

The above stated formula can be decomposed and expressed also as

− −

{ }

1

−

−

=−∞ =−∞ =−∞

In the papers (Brančík, 2007a, 2007b, 2010b), it was shown for the 1D, 2D, and 3D cases, the rectangular method of a numerical integration leads to an approximate formula whose a relative error is adjustable, and corresponds to the complex Fourier series approximation of a respective dimension. The method has been generalized for an arbitrary dimension *n* in the recent work (Brančík, 2011).

#### **2.1 Complex Fourier series approximation and limiting relative error**

Substituting *si* = *ci* + *jωi* into (2), and using a rectangular rule of the integration, namely ω*i* = *mi*Ω*i*, and Ω*i* = 2π/τ*i* as generalized frequency steps, with τ*i* forming a region of the solution *t* ∈ [0,τ1) ×...× [0,τ*n*), an approximate formula is

$$\tilde{f}(\mathbf{t}) = \exp(\mathbf{c}\mathbf{t}^T) \left(\prod\_{i=1}^n \mathbf{r}\_i^{-1}\right) \sum\_{m\_1 = -\infty}^\infty \cdots \sum\_{m\_n = -\infty}^\infty F(\mathbf{s}) \exp\left(j \sum\_{i=1}^n m\_i \Omega\_i t\_i\right) \tag{3}$$

with *si* = *ci* + *jmi*Ω*i*, ∀*i* . As is shown in (Brančík, 2011), a limiting relative error δ*M* of (3) can be controlled by setting **c** = (*c*1,...,*cn*), defining paths of the integration in (2), namely

$$c\_i = \alpha\_i - \frac{1}{\pi\_i} \ln \left( 1 - \frac{1}{\sqrt[n]{1 + \delta\_M}} \right) = \alpha\_i - \frac{1}{\pi\_i} \ln \frac{\delta\_M}{n} \,\tag{4}$$

for *i* = 1,…,*n*, and while keeping the equalities τ1(*c*1 – α1) =...= τ*n*(*cn* – α*n*). The simplification in (4) is enabled due to small values δ*M* considered in practice. The last equation is used for setting up parameters of the *n*D NILT method relating them to a limiting relative error δ*M* required for practical computations.

#### **2.2 Practical computational methods**

It should be highlighted that the formula (4) is valid, and a relative error supposed is really achievable by the *n*D NILT (3), if infinite numbers of terms are used in the series. In practice, it cannot sure be fulfilled, but a suitable technique for accelerating a convergence of infinite series is usable, as is e.g. a quotient-difference (q-d) algorithm (Rutishauser, 1957). Besides, as has been already successfully used for cases of *n* ≤ 3, the formula (3) can be rearranged to enable using FFT & IFFT algorithms for an effective computation.

#### **2.2.1 Partial ILTs evaluation technique**

The technique of practical evaluation of the *n*-fold infinite sum (3) follows the properties of the *n*-fold Bromwich integral (2), namely we can rearrange it into the form

$$f(t\_1, t\_2, \ldots, t\_n) = \frac{1}{2\pi j} \int\_{c\_1 - j\alpha}^{c\_1 + j\alpha} \left| \frac{1}{2\pi j} \int\_{c\_2 - j\alpha}^{c\_2 + j\alpha} \right|^{\alpha} \cdots \frac{1}{2\pi j} \int\_{c\_n - j\alpha}^{c\_n + j\alpha} F(s\_1, s\_2, \ldots, s\_n) e^{s\_n t\_n} ds\_n \cdots \right| e^{s\_1 t\_2} ds\_2 \Big| e^{s\_1 t\_1} ds\_1 \wedge \ldots \wedge e^{s\_n t\_n} ds\_n \wedge \Big|$$

or shortly

$$f(\mathbf{t}\_1, \mathbf{t}\_2, \dots, \mathbf{t}\_n) = \mathbb{L}\_1^{-1} \left[ \mathbb{L}\_2^{-1} \left[ \cdots \mathbb{L}\_n^{-1} [F(\mathbf{s}\_1, \mathbf{s}\_2, \dots, \mathbf{s}\_n)] \cdots \right] \right]. \tag{6}$$

Although the order of the integration may be arbitrary on principle, here the above one will be used for an explanation. Similarly, (3) can be rewritten as

In the papers (Brančík, 2007a, 2007b, 2010b), it was shown for the 1D, 2D, and 3D cases, the rectangular method of a numerical integration leads to an approximate formula whose a relative error is adjustable, and corresponds to the complex Fourier series approximation of a respective dimension. The method has been generalized for an arbitrary dimension *n* in

Substituting *si* = *ci* + *jωi* into (2), and using a rectangular rule of the integration, namely

1 1

*i i*

 α

> τ

*f t ct s j mt* , (3)

 Ω

*n*

*i i=*

*n n*

τ

*M*

δ

*n*

*M* considered in practice. The last equation is used for

*n n*

*i i*

*i* forming a region of the

δ

, (4)

2 2 1 1

*s t st st*

*n*). The simplification

*M* of (3) can

δ*M*

*i* as generalized frequency steps, with

*n*), an approximate formula is

1

with *si* = *ci* + *jmi*Ω*i*, ∀*i* . As is shown in (Brančík, 2011), a limiting relative error

be controlled by setting **c** = (*c*1,...,*cn*), defining paths of the integration in (2), namely

τ

( ) = exp( ) ( )exp

1

∏

*m m*

11 1 = ln 1 ln 1

δ

τ1(*c*1 – α1) =...= τ*n*(*cn* – α

 − − ≈− +

setting up parameters of the *n*D NILT method relating them to a limiting relative error

It should be highlighted that the formula (4) is valid, and a relative error supposed is really achievable by the *n*D NILT (3), if infinite numbers of terms are used in the series. In practice, it cannot sure be fulfilled, but a suitable technique for accelerating a convergence of infinite series is usable, as is e.g. a quotient-difference (q-d) algorithm (Rutishauser, 1957). Besides, as has been already successfully used for cases of *n* ≤ 3, the formula (3) can be rearranged to

The technique of practical evaluation of the *n*-fold infinite sum (3) follows the properties of

1 2 1 2 2 1

 , (5)

1 2 1 2 1 2 ( , , , )= (, , , ) *n nn ft t t Fs s s* −− − . (6)

[ ] 11 1

*n*

*f tt t F s s s e ds e ds e ds*

*n*

Although the order of the integration may be arbitrary on principle, here the above one will

*n n n*

 π

*i i i n i i M*

 *F* ∞ ∞ − = =−∞ =−∞

**2.1 Complex Fourier series approximation and limiting relative error**

the recent work (Brančík, 2011).

τ

π/τ

τ

*c*

for *i* = 1,…,*n*, and while keeping the equalities

in (4) is enabled due to small values

required for practical computations.

**2.2 Practical computational methods**

**2.2.1 Partial ILTs evaluation technique**

or shortly

α

τ

enable using FFT & IFFT algorithms for an effective computation.

1 2

11 1 22 2 ( , , , )= (, , , )

*jj j*

1 2

be used for an explanation. Similarly, (3) can be rewritten as

ππ

the *n*-fold Bromwich integral (2), namely we can rearrange it into the form

+∞ +∞ + ∞

−∞ −∞ − ∞

*cj c j c j*

*cj c j c j*

δ

*T*

1) ×...× [0,

*i* = *mi*Ω*i*, and Ω*i* = 2

solution *t* ∈ [0,

ω

$$\tilde{f}\left(t\_1, t\_2, \ldots, t\_n\right) = \frac{e^{\varepsilon\_1 t\_1}}{\tau\_1} \sum\_{m\_1 = -\alpha}^{\alpha} \left(\frac{e^{\varepsilon\_2 t\_2}}{\tau\_2} \sum\_{m\_2 = -\alpha}^{\alpha} \left(\cdots \frac{e^{\varepsilon\_n t\_n}}{\tau\_n} \sum\_{m\_n = -\alpha}^{\alpha} F(s\_1, s\_2, \ldots, s\_n) e^{jm\_n \Omega\_n t\_n} \cdots \right) e^{im\_2 \Omega\_2 t\_2} \right) e^{im\_1 \Omega\_1 t\_1} \ , \ \{\mathcal{T}\}$$

with *si* = *ci* + *jmi*Ω*i*. If we define *Fn* ≡ *F*(*s*1,...,*sn*-1,*sn*) and *F*0 ≡ *f*(*t*1,...,*tn*-1,*tn*), then *n* consequential partial inversions are performed as

$$\begin{aligned} \mathbb{L}\_n^{-1} \{ F\_n \} &= F\_{n-1} (\mathbf{s}\_1, \dots, \mathbf{s}\_{n-1}, t\_n) \\ \mathbb{L}\_{n-1}^{-1} \{ F\_{n-1} \} &= F\_{n-2} (\mathbf{s}\_1, \dots, \mathbf{t}\_{n-1}, t\_n) \\ &\vdots \\ \mathbb{L}\_1^{-1} \{ F\_1 \} &= f(t\_1, \dots, t\_{n-1}, t\_n) \end{aligned} \tag{8}$$

As is obvious we need to use a procedure able to make the inversion of Laplace transforms dependent on another *n-*1 parameters, complex in general. Let us denote arguments in (8) by *pi* = (*p*1,...,*pn*-1,*pn*). Then the ILT of the type

$$F\_{i-1}(\mathcal{p}\_{i-1}) = \mathbb{L}\_i^{-1} \left\{ F\_i(\mathcal{p}\_i) \right\} = \frac{1}{2\pi j} \int\_{c\_i - j^{\rm iso}}^{c\_i + j^{\rm iso}} F\_i(\mathcal{p}\_i) e^{s\_i t\_i} \, ds\_i \tag{9}$$

can be used *n* times, *i* = *n*,*n*-1...,1, to evaluate (8), with *pn* = (*s*1,...,*sn*-1,*sn*), *pn*-1 = (*s*1,...,*sn*-1,*tn*),..., *p*1 = (*s*1,...,*tn*-1,*tn*), and *p*0 = (*t*1,...,*tn*-1,*tn*), while *pj* = *sj* for *j* ≤ *i*, and *pj* = *tj* otherwise. A further technique is based on demand to find the solution on a whole region of discrete points. Then, taking into account *tik* = *kTi* in (9), with *Ti* as the sampling periods in the original domain, we can write an approximate formula

$$\tilde{F}\_{i-1}(\mathcal{p}\_{i-1}) = \frac{e^{c\_i kT\_i}}{\pi\_i} \sum\_{m = -\infty}^{\infty} \tilde{F}\_i(\mathcal{p}\_i) e^{j2\pi mkT\_i f \mathbf{r}\_i} \quad , \tag{10}$$

*i* = *n*,*n*-1...,1, and with Ω*i* = 2π/τ*i* substituted. As follows from the error analysis (Brančík, 2011) a relative error is predictable on the region Οerr = [0,τ1) ×...× [0,τ*n*). For *k* = 0*,*1*,...,Mi*-1, *i* = 1,...,*n*, a maximum reachable region is Οmax = [0,(*M*1-1)*T*1] ×...× [0,(*Mn*-1)*Tn*]. Thus, to meet the necessary condition Οmax ⊂ Οerr, we can set up fittingly τ*i* = *MiTi*, *i* = 1,...,*n*. In practice, a region of the calculation is chosen to be Οcal = [0,*t*1cal] ×...× [0,*tn*cal], with *ti*cal = (*Mi*/2–1)*Ti*, *i* = 1,...,*n*, to provide certain margins.

#### **2.2.2 FFT, IFFT, and quotient-difference algorithms utilization**

As is shown in (Brančík, 2007a, 2010c), the discretized formula (10) can be evaluated by the FFT and IFFT algorithms, in conjunction with the quotient-difference (q-d) algorithm for accelerating convergence of the residual infinite series, see following procedures.

To explain it in more detail, let us consider an *r*-th cycle in gaining the original function via (9), i.e. { } 1 1 1 ( )= ( ) *F F r r r rr* − − − *p p* . For its discretized version (10) we have

$$\tilde{F}\_{r-1}(\mathbf{s}\_1, \dots, k T\_r, \dots, t\_n) = \frac{e^{c\_r k \tilde{T}\_r}}{\pi\_r} \sum\_{m = -\infty}^{\infty} \tilde{F}\_r(\mathbf{s}\_1, \dots, \mathbf{c}\_r + j m \frac{2\pi}{\tau\_r}, \dots, t\_n) e^{j2\pi m k T\_r f \tau\_r} \ . \tag{11}$$

The above stated formula can be decomposed and expressed also as

$$\tilde{F}\_{r-1}(\mathbf{s}\_1, \dots, k\boldsymbol{\Gamma}\_r, \dots, \mathbf{t}\_n) = \frac{e^{c\_r kT\_r}}{\pi\_r} \left[ \sum\_{m=0}^{M\_r - 1} \tilde{F}\_r^{(-m)} \mathbf{z}\_{-k}^m + \sum\_{m=0}^{\alpha} \tilde{\mathbf{G}}\_r^{(-m)} \mathbf{z}\_{-k}^m + \sum\_{m=0}^{M\_r - 1} \tilde{F}\_r^{(m)} \mathbf{z}\_k^m + \sum\_{m=0}^{\alpha} \tilde{\mathbf{G}}\_r^{(m)} \mathbf{z}\_k^m - \tilde{\mathbf{F}}\_r^{(0)} \right],\\\text{(12)}$$

where individual terms are defined as

$$\begin{aligned} M\_r &= \mathbb{Z}^{k\_r}, K\_r \text{ integer }, \\ \tilde{F}\_r^{(\pm m)} &= \tilde{F}\_r(s\_1, \dots, c\_r \pm jm2\pi/\mathfrak{r}\_r, \dots, t\_n) \\ \tilde{G}\_r^{(\pm m)} &= \tilde{F}\_r^{(\pm M\_r \pm m)}, \\ z\_{\pm k} &= \exp(\pm j2\pi kT\_r/\mathfrak{r}\_r) \end{aligned} \tag{13}$$

Numerical Inverse Laplace Transforms for Electrical Engineering Simulation 55

For practical computations, however, the recursive formulae stated below are more effective

for *m =* 1*,...,2P*, ∀*k* , with the initial values *A*-1 = 0, *B*-1 = 1, *A*0 = *d*0, and *B*0 = 1. Then, instead of

2 2

*G z Az Bz k*

± ±±

( ) ( ) ( ), *m m r k Pk Pk*

The q-d algorithm is a very efficient tool just for a power series convergence acceleration, here enabling (7) to achieve a relative error near its theoretical value defined by (4), see the

In this part experimental verifications of the *n*D NILT theory above will first be presented, for one to three dimensional cases, i.e. *n* ≤ 3. For such dimensions the Matlab functions have been developed and errors stated on a basis of some sample images with known originals. The Matlab listings of basic versions of the NILT functions are provided, together with examples of their right calling. Another Matlab listings will be discussed in more detail later,

In case of the 1D inverse LT, a well-known Bromwich integral results from (2), namely

1 ( )= ( ) 2 *c st f t F s e dt* + ∞

− ∞ π

where indexes 1 were omitted. By using the theory above a path of the numerical

1 1 11 1 1 1

= ln = ln ln *M M M*

In contrast to most other approaches, the 1D NILT method described here enables to treat complex images resulting in complex originals as no real or imaginary parts are extracted during an evaluation process. It can be useful in some special applications, not only in the

1 = = ( ) = = cos sin *s j t F s j f te t j t*

+ + − + +

ω

> ω

τδ

α

 − − + + ≈− +

1

electrical engineering. We can śhow it on a simple transform pair

ω

( ) 22 22

*s j s s*

ω  δ *j*

1 2 12 *A z A z dz A z B z B z dz B z m k m k m k m k m k m k m km k* ( )= ( ) ( ), ( )= ( ) ( ) ± −± ± − ± ± −± ± − ± + + , (18)

0 2 1 2 1 ( ) ( ) ( ) = , = , = , =, , *rj j j j dG d* − − − *q d e j P* . (17)

≈ ∀ . (19)

*c j j* , (20)

 αδ

ω  τ

ω

. (22)

 ω

. (21)

0 0 0

Then, the coefficients *dm*, *m* = 0,...,2*P*, in (14) are given by

to be used (DeHoog et al., 1982). They are of the forms

0

=

**2.3 Matlab listings and experimental errors evaluation**

∞ ±

*m*

the continued fraction (14), we can write

in the chapter with practical applications.

integration is stated according to (4), leading to

τ

*c* α

**2.3.1 One-dimensional NILT**

following examples.

when 2 = =, 1 *Mr j k k ze k* ± π ± ∀ , has been considered, and τ*r* = *MrTr*. As is evident the first and the third finite sum of (12) can be evaluated via the FFT and IFFT algorithms, respectively, while 2*P*+1 terms from the infinite sums are used as the input data in the quotient-difference algorithm (Macdonald, 1964; McCabe, 1983; Rutishauser, 1957). We can replace the above infinite power series by a continued fraction as

$$\sum\_{m=0}^{\infty} \tilde{G}\_r^{(\pm m)} z\_{\pm k}^m = d\_0 \left\langle 1 + d\_1 z\_{\pm k} \left\langle 1 + \dots + d\_{2^p} z\_{\pm k} \right\rangle \right\rangle \,' \forall k \,\,'\tag{14}$$

which gives much more accurate result than the original sum truncated on 2*P*+1 terms only. The q-d algorithm process can be explained based on a lozenge diagram shown in Fig. 1.

Fig. 1. Quotient-difference algorithm lozenge diagram

The first two columns are formed as

$$\begin{aligned} e\_0^{(i)} &= 0 \\ q\_1^{(i)} &= \tilde{G}\_r^{\pm(i+1)} \Big/ \tilde{G}\_r^{\pm i} \quad \text{i} = 0, \ldots, 2P - 1 \end{aligned} \tag{15}$$

while the successive columns are given by the rules

$$\begin{aligned} e\_j^{(i)} &= q\_j^{(i+1)} - q\_j^{(i)} + e\_{j-1}^{(i+1)}, \quad i = 0, \dots, 2P - 2j \ &\quad \text{for} \quad j = 1, \dots, P \ &\quad \text{(16)}\\ q\_j^{(i)} &= q\_{j-1}^{(i+1)} e\_{j-1}^{(i+1)} \int e\_{j-1}^{(i)} \ &\quad \text{i = 0, \dots, 2P - 2j - 1, \quad \text{for} \quad j = 2, \dots, P \ &\quad \end{aligned} \tag{16}$$

Then, the coefficients *dm*, *m* = 0,...,2*P*, in (14) are given by

54 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

1 1

, (12)

*r m m mm*

( ) ( ) ( ) ( ) () ( , , , , )= *c kT r r M M r r*

*r rn r k r k r k rk r*

*e F s kT t F z G z Fz Gz F*

1

= ,

*z j kT*

*r*

2

*k r r*

±

= exp( ) ,

π

= , integer ,

2

*M K*

*r*

*K r r*

( ) ( )

± ± ±

We can replace the above infinite power series by a continued fraction as

0 0

( )

*e*

4 0

( )

*e*

( ) ( )

1

−

Fig. 1. Quotient-difference algorithm lozenge diagram

while the successive columns are given by the rules

1 1 11 1

+ + −− −

() ( ) () ( ) () ( ) ( ) ()

*i i ii j j jj i ii i j jj j*

0

( )

*i*

1

1 1

+ +

The first two columns are formed as

*m M m r r*

τ

where individual terms are defined as

− − −

( )

±

*m*

*G F*

±

± ∀ , has been considered, and

0

=

∞ ±

*m*

0 0 00

= = ==

2

π τ

τ

= ( , , , , ),

 τ

*r r r rn*

±

As is evident the first and the third finite sum of (12) can be evaluated via the FFT and IFFT algorithms, respectively, while 2*P*+1 terms from the infinite sums are used as the input data in the quotient-difference algorithm (Macdonald, 1964; McCabe, 1983; Rutishauser, 1957).

01 2

1 1 ( ) ( ( )) *m m r k k P k*

which gives much more accurate result than the original sum truncated on 2*P*+1 terms only. The q-d algorithm process can be explained based on a lozenge diagram shown in Fig. 1.

> 1 0 1 2 210 012 2 1 1 2

*q q e ee q q*

( ) ( ) ( ) ( ) ( ) ( ) ( )

±± ±

*G z d dz d z*

0 1 1 0 0 1

( ) ( ) ( )

*q e e*

3 2 0 1 3 1

() () ( )

*e e q*

1

− + −

*i i i r r*

0 0 2

*e iP qG G i P* ± ± + −

= , =, , ,

= , =, , ,

= , = , , , for = , , , = , = , , , for = , , .

*e q q e i Pj j P q qe e i Pj j P*

0 21

022 1 0 2 21 2

(16)

− −

(15)

*F F s c jm t*

− − ∞ ∞ − −

*m m mm m m mm*

(13)

*r* = *MrTr*.

≈ + ++ , ∀*k* , (14)

+ + +−

1 1

when 2

= =, 1 *Mr j k k ze k* ± π

0

$$d\_0 = \tilde{G}\_r^{(0)} \; , \; d\_{2j-1} = -q\_j^{(0)} \; , \; \; d\_{2j} = -e\_j^{(0)} \; , \; \; j = 1 \; \; \dots \; \; P \; . \tag{17}$$

For practical computations, however, the recursive formulae stated below are more effective to be used (DeHoog et al., 1982). They are of the forms

$$A\_m(z\_{\pm k}) = A\_{m-1}(z\_{\pm k}) + d\_m z\_{\pm k} A\_{m-2}(z\_{\pm k}) \ , \ B\_m(z\_{\pm k}) = B\_{m-1}(z\_{\pm k}) + d\_m z\_{\pm k} B\_{m-2}(z\_{\pm k}) \ , \tag{18}$$

for *m =* 1*,...,2P*, ∀*k* , with the initial values *A*-1 = 0, *B*-1 = 1, *A*0 = *d*0, and *B*0 = 1. Then, instead of the continued fraction (14), we can write

$$\sum\_{m=0}^{\infty} \tilde{\mathbf{G}}\_r^{(\pm m)} z\_{\pm k}^m = A\_{2P}(z\_{\pm k}) / \mathcal{B}\_{2P}(z\_{\pm k}) \; , \; \forall k \; . \tag{19}$$

The q-d algorithm is a very efficient tool just for a power series convergence acceleration, here enabling (7) to achieve a relative error near its theoretical value defined by (4), see the following examples.

#### **2.3 Matlab listings and experimental errors evaluation**

In this part experimental verifications of the *n*D NILT theory above will first be presented, for one to three dimensional cases, i.e. *n* ≤ 3. For such dimensions the Matlab functions have been developed and errors stated on a basis of some sample images with known originals. The Matlab listings of basic versions of the NILT functions are provided, together with examples of their right calling. Another Matlab listings will be discussed in more detail later, in the chapter with practical applications.

#### **2.3.1 One-dimensional NILT**

In case of the 1D inverse LT, a well-known Bromwich integral results from (2), namely

$$f(t) = \frac{1}{2\pi j} \int\_{c - j^{\rm os}}^{c + j^{\rm iso}} F(s)e^{st} \, dt \,\,\,\,\,\tag{20}$$

where indexes 1 were omitted. By using the theory above a path of the numerical integration is stated according to (4), leading to

$$\alpha = \alpha - \frac{1}{\pi} \ln \left( 1 - \frac{1}{1 + \mathcal{S}\_{\mathcal{M}}} \right) = \alpha + \frac{1}{\pi} \ln \left( 1 + \frac{1}{\mathcal{S}\_{\mathcal{M}}} \right) = \alpha - \frac{1}{\pi} \ln \mathcal{S}\_{\mathcal{M}} \,. \tag{21}$$

In contrast to most other approaches, the 1D NILT method described here enables to treat complex images resulting in complex originals as no real or imaginary parts are extracted during an evaluation process. It can be useful in some special applications, not only in the electrical engineering. We can śhow it on a simple transform pair

$$F(s) = \frac{l}{s - j\alpha} = \frac{s}{s^2 + \alpha^2} + j\frac{\alpha\alpha}{s^2 + \alpha^2} \quad \mapsto \quad f(t) = e^{i\alpha t} = \cos\alpha t + j\sin\alpha t \,\,. \tag{22}$$

Numerical Inverse Laplace Transforms for Electrical Engineering Simulation 57

Graphical results and corresponding errors are shown in Fig. 2. Because the originals are

Im[f(t)]

= 0, we can see the errors satisfy (21) very well (

10-15

10-10

10-5

Error

Another test functions are considered in Tab. 2, with numerical results shown in Fig. 3. As is again obvious from Fig. 3 the relative errors satisfy theoretical expectations, with an

i 1 2 3 4 5 6

π

2

1 *s* +1

*t* 0*J t*( ) 1

11 2 2

δ

− , (23)

*st s t*

+

2 *s* 4 π

+ π

In case of the 2D inverse LT, a two-fold Bromwich integral results from (2), namely

1 2

+∞ +∞

−∞ −∞

and by using the theory above the paths of numerical integrations are stated based on (4) as

1 4 ( , )= (, ) *c jc j*

π

*i i i*

1 2

*c jc j f t t F s s e ds ds*

*i i M c i*

 − − ≈− +

δ

1 2 2 1 2 1 2

11 1 1 1 2 1 2 = ln ln , = , *M*

 τ

 α

100

0 1 2 3 4

Imaginary part of the original

t

0 1 2 3 4

t

*s e s*

2 erfc

*t*

− *s e*

1( ) *t* −1

*s* −

. (24)

δ

*M* = 10-10 was

bounded by values ±1, and

10-15

exception of vicinities of discontinuities.

*s* +1 2

Table 2. Test Laplace transforms for errors evaluation

α

τ

1

10-10

10-5

Error

( ) *F s i*

( ) *i f t t e*

**2.3.2 Two-dimensional NILT**

100

Re[f(t)]

α

0 1 2 3 4

Real part of the original

t

0 1 2 3 4

t

1

− *t te*− sin( ) 2

Fig. 2. Numerical inversion leading to complex original *f*(*t*) = exp(*jωt*)

( ) *s* +1 2 2

considered), excluding only beginning of the interval.

Of course, when preprocessing the transform to arrange it to a Cartesian form, as is shown on the right sides in (22), the result could by get by inverting the real and imaginary parts separately, by using an arbitrary NILT method. Here, however, no symbolic manipulations are needed in advance, and *F*(*s*) enters the NILT function in its basic form as a whole.

A Matlab language listing is shown in Tab. 1, where the relative error needed is marked by Er and is subject to a change if necessary, similarly as the minimal abscissa of convergence (exponential order) α , alfa, numbers of points for the resultant solution, M, and for the q-d algorithm, P. If only real transforms *F*(*s*) are considered the bottom line in the listing can be inactivated. The NILT function is called from a command line as follows: niltc('F',tm); where F is a name of another function in which the *F*(*s*) is defined, and tm marks an upper limit of the original variable *t*. In our case, and for ω = 2π, this function can have a form

```
function f=expc(s)
f=1./(s-2*pi*j);
```

```
% ------ 1D NILT for complex arguments – basic version -------
% ----------- based on FFT/IFFT/q-d, by L. Brančík -----------
function [ft,t]=niltc(F,tm);
alfa=0; M=256; P=3; Er=1e-10; % adjustable
N=2*M; qd=2*P+1;
t=linspace(0,tm,M); NT=2*tm*N/(N-2); omega=2*pi/NT;
c=alfa+log(1+1/Er)/NT; s=c-i*omega*(0:N+qd-1);
Fs(1,:)=feval(F,s); Fs(2,:)=feval(F,conj(s));
ft(1,:)=fft(Fs(1,1:N)); ft(2,:)=N*ifft(Fs(2,1:N));
ft=ft(:,1:M); D=zeros(2,qd); E=D;
Q=Fs(:,N+2:N+qd)./Fs(:,N+1:N+qd-1);
D(:,1)=Fs(:,N+1); D(:,2)=-Q(:,1);
for r=2:2:qd-1
 w=qd-r;
 E(:,1:w)=Q(:,2:w+1)-Q(:,1:w)+E(:,2:w+1); D(:,r+1)=-E(:,1);
 if r>2
 Q(:,1:w-1)=Q(:,2:w).*E(:,2:w)./E(:,1:w-1);
       D(:,r)=-Q(:,1);
 end
end
A2=zeros(2,M); B2=ones(2,M); A1=repmat(D(:,1),[1,M]); B1=B2;
z1=exp(-i*omega*t); z=[z1;conj(z1)];
for n=2:qd
 Dn=repmat(D(:,n),[1,M]);
 A=A1+Dn.*z.*A2; B=B1+Dn.*z.*B2; A2=A1; B2=B1; A1=A; B1=B;
end
ft=ft+A./B; ft=sum(ft)-Fs(2,1); ft=exp(c*t)/NT.*ft;
ft(1)=2*ft(1);
figure; plot(t,real(ft));
figure; plot(t,imag(ft)); % optional
```
Table 1. Matlab listing of 1D NILT method accepting complex arguments

As is obvious, the Laplace transform must be defined to enable Matlab array processing, i.e. element-by-element array operators have to be used. Thus, the calling our function can look like niltc('expc',4); if the function is saved under the same name, expc, or it is placed inside the M-file with own NILT function (Tab. 1), following always its body. Alternatively, the calling can look like [ft,t]=niltc('F',tm); if respective variables in the brackets are to be saved in the memory after the function finishes.

Of course, when preprocessing the transform to arrange it to a Cartesian form, as is shown on the right sides in (22), the result could by get by inverting the real and imaginary parts separately, by using an arbitrary NILT method. Here, however, no symbolic manipulations

A Matlab language listing is shown in Tab. 1, where the relative error needed is marked by Er and is subject to a change if necessary, similarly as the minimal abscissa of convergence

algorithm, P. If only real transforms *F*(*s*) are considered the bottom line in the listing can be inactivated. The NILT function is called from a command line as follows: niltc('F',tm); where F is a name of another function in which the *F*(*s*) is defined, and tm marks an upper

% ------ 1D NILT for complex arguments – basic version ------- % ----------- based on FFT/IFFT/q-d, by L. Brančík -----------

alfa=0; M=256; P=3; Er=1e-10; % adjustable

E(:,1:w)=Q(:,2:w+1)-Q(:,1:w)+E(:,2:w+1); D(:,r+1)=-E(:,1);

A2=zeros(2,M); B2=ones(2,M); A1=repmat(D(:,1),[1,M]); B1=B2;

A=A1+Dn.\*z.\*A2; B=B1+Dn.\*z.\*B2; A2=A1; B2=B1; A1=A; B1=B;

figure; plot(t,imag(ft)); % optional

t=linspace(0,tm,M); NT=2\*tm\*N/(N-2); omega=2\*pi/NT; c=alfa+log(1+1/Er)/NT; s=c-i\*omega\*(0:N+qd-1); Fs(1,:)=feval(F,s); Fs(2,:)=feval(F,conj(s)); ft(1,:)=fft(Fs(1,1:N)); ft(2,:)=N\*ifft(Fs(2,1:N));

Q(:,1:w-1)=Q(:,2:w).\*E(:,2:w)./E(:,1:w-1);

ft=ft+A./B; ft=sum(ft)-Fs(2,1); ft=exp(c\*t)/NT.\*ft;

ω = 2π

, alfa, numbers of points for the resultant solution, M, and for the q-d

, this function can have a form

are needed in advance, and *F*(*s*) enters the NILT function in its basic form as a whole.

Table 1. Matlab listing of 1D NILT method accepting complex arguments

are to be saved in the memory after the function finishes.

As is obvious, the Laplace transform must be defined to enable Matlab array processing, i.e. element-by-element array operators have to be used. Thus, the calling our function can look like niltc('expc',4); if the function is saved under the same name, expc, or it is placed inside the M-file with own NILT function (Tab. 1), following always its body. Alternatively, the calling can look like [ft,t]=niltc('F',tm); if respective variables in the brackets

(exponential order)

function f=expc(s) f=1./(s-2\*pi\*j);

α

limit of the original variable *t*. In our case, and for

function [ft,t]=niltc(F,tm);

ft=ft(:,1:M); D=zeros(2,qd); E=D; Q=Fs(:,N+2:N+qd)./Fs(:,N+1:N+qd-1); D(:,1)=Fs(:,N+1); D(:,2)=-Q(:,1);

D(:,r)=-Q(:,1);

Dn=repmat(D(:,n),[1,M]);

figure; plot(t,real(ft));

z1=exp(-i\*omega\*t); z=[z1;conj(z1)];

N=2\*M; qd=2\*P+1;

for r=2:2:qd-1 w=qd-r;

if r>2

for n=2:qd

ft(1)=2\*ft(1);

 end end

end

Graphical results and corresponding errors are shown in Fig. 2. Because the originals are bounded by values ±1, and α = 0, we can see the errors satisfy (21) very well (δ*M* = 10-10 was considered), excluding only beginning of the interval.

Fig. 2. Numerical inversion leading to complex original *f*(*t*) = exp(*jωt*)

Another test functions are considered in Tab. 2, with numerical results shown in Fig. 3. As is again obvious from Fig. 3 the relative errors satisfy theoretical expectations, with an exception of vicinities of discontinuities.

Table 2. Test Laplace transforms for errors evaluation

#### **2.3.2 Two-dimensional NILT**

In case of the 2D inverse LT, a two-fold Bromwich integral results from (2), namely

$$f(t\_1, t\_2) = -\frac{1}{4\pi^2} \int\_{c\_1 - j\alpha}^{c\_1 + j\alpha} \int\_{-j\alpha}^{j\alpha} F(s\_1, s\_2) e^{s\_1 t\_1 + s\_2 t\_2} ds\_1 ds\_2 \,\tag{23}$$

and by using the theory above the paths of numerical integrations are stated based on (4) as

$$\alpha c\_i = \alpha\_i - \frac{1}{\pi\_i} \ln \left( 1 - \frac{1}{\sqrt{1 + \delta\_M}} \right) = \alpha\_i - \frac{1}{\pi\_i} \ln \frac{\delta\_M}{2}, \quad i = 1, 2 \dots \tag{24}$$

Numerical Inverse Laplace Transforms for Electrical Engineering Simulation 59

% ----- 2D NILT based on partial inversions, by L. Brančík -----

alfax=0; alfay=0; Mx=256; My=256; P=3; Er=1e-8; % adjustable xpl=64; ypl=64; % adjustable

x=linspace(0,xm,Mx); y=linspace(0,ym,My); x=x(rx); y=y(ry); [q,p]=meshgrid(Asigy,Asigx2); Fpq(:,:,1)=feval(F,p,q); [q,p]=meshgrid(conj(Asigy),Asigx2); Fpq(:,:,2)=feval(F,p,q); Fpyp=Pnilt(Fpq,Ny,ry,qd,y,ny,omegay,sigy); % Pnilt to get F(p,y)

fxy=Pnilt(Fpy,Nx,rx,qd,x,nx,omegax,sigx); % Pnilt to get f(x,y)

figure; mesh(x,y,imag(fxy)); % optional % ------ PARTIAL NILT based on FFT/IFFT/Q-D, by L.Brančík ------

fx(:,:,1)=fft(Fq(:,:,1),N,2); fx(:,:,2)=N\*ifft(Fq(:,:,2),N,2);

w=qd-r; e(:,1:w,:)=q(:,2:w+1,:)-q(:,1:w,:)+e(:,2:w+1,:);

q(:,1:w-1,:)=q(:,2:w,:).\*e(:,2:w,:)./e(:,1:w-1,:);

z1(1,:,1)=exp(-i\*omega\*xy); z1(1,:,2)=conj(z1(1,:,1));

fx=fx+A./B; fx=sum(fx,3)-repmat(Fq(:,1),[1,delxy]); fx=repmat(exp(c\*xy)/nxy,[delv,1]).\*fx; fx(:,1)=2\*fx(:,1);

d=zeros(delv,qd,2); e=d; q=Fq(:,N+2:N+qd,:)./Fq(:,N+1:N+qd-1,:);

Nx=2\*Mx; Ny=2\*My; qd=2\*P+1; Ke=log(1-1/sqrt(1+Er));

omegax=2\*pi/nx; omegay=2\*pi/ny; sigx=alfax-Ke/nx; sigy=alfay-Ke/ny; qd1=qd-1; Nxw=Nx+qd1; Nyw=Ny+qd1; Asigx=sigx-i\*omegax\*(0:Nxw); Asigy=sigy-i\*omegay\*(0:Nyw);

function fx=Pnilt(Fq,N,grid,qd,xy,nxy,omega,c);

A2=zeros(delv,delxy,2); B2=ones(delv,delxy,2);

A1=repmat(d(:,1,:),[1,delxy,1]); B1=B2;

Dn=repmat(d(:,n,:),[1,delxy,1]);

d(:,1,:)=Fq(:,N+1,:); d(:,2,:)=-q(:,1,:);

fx=fx(:,grid,:); delv=size(Fq,1); delxy=length(xy);

function fxy=nilt2c(F,xm,ym);

nx=2\*xm\*Nx/(Nx-2); ny=2\*ym\*Ny/(Ny-2);

rx=[1:Mx/xpl:Mx,Mx]; ry=[1:My/ypl:My,My];

Asigx2=cat(2,Asigx,conj(Asigx));

Fpy(:,:,1)=Fpyp(1:Nxw+1,:).'; Fpy(:,:,2)=Fpyp(Nxw+2:2\*Nxw+2,:).';

figure; mesh(x,y,real(fxy));

d(:,r+1,:)=-e(:,1,:);

d(:,r,:)=-q(:,1,:);

z=repmat(z1,[delv,1]);

Table 3. Matlab listing of 2D NILT based on partial inversions

( ) 1 2

absolute error. The respective Matlab function can be of a form

1 2 exp( ) 2

*s s*

corresponding transform pair is

for r=2:2:qd-1

if r>2

for n=2:qd

end

 end end

Another simple example shows a shifted 2D unit step, with different shifts along the axis. A

A=A1+Dn.\*z.\*A2; B=B1+Dn.\*z.\*B2; A2=A1; B2=B1; A1=A; B1=B;

1 2 12 1 2

, = ( , ) = 1( 2, 1) *s s Fs s ft t t t*

In this case, a displaying imaginary part gives a zero plane, and the respective line in the 2D NILT function can be inactivated. The graphical results are depicted in Fig. 5, including an

− − − − . (27)

Fig. 3. Computed originals and errors for test Laplace transforms in Tab. 2

A Matlab language listing is shown in Tab. 3, with all the parameters denoted by similar way as in the previous case.

Nevertheless, the Laplace transform variables and the original variables have changed in this listing as *s*1 → *p*, *s*2 → *q*, and *t*1 → *x*, *t*2 → *y*, respectively, which simplified a writing. The parameters are then indexed in compliance with these new notations. Besides, numbers of points used to plot three-dimensional graphical results are set by xpl and ypl.

For the same reasons as at the 1D NILT, the 2D NILT method discussed here enables to treat complex images of two variables resulting in complex originals. We will śhow it on a simple transform pair

$$F(s\_1, s\_2) = \frac{1}{(s\_1 - jo\_1)(s\_2 - jo\_2)} \quad \mapsto \quad f(t\_1, t\_2) = e^{j(o\_1t\_1 + o\_2t\_2)}.\tag{25}$$

After rearranging the above equation, we can also write

$$F(s\_1, s\_2) = \frac{s\_1 s\_2 - a \rho\_1 \rho\_2}{\left(s\_1^2 + a \rho\_1^2\right)\left(s\_2^2 + a \rho\_2^2\right)} + j \frac{a \rho\_2 s\_1 + a \rho\_1 s\_2}{\left(s\_1^2 + a \rho\_1^2\right)\left(s\_2^2 + a \rho\_2^2\right)} \quad \leftrightarrow$$

$$f(t\_1, t\_2) = \cos\left(a\rho\_1 t\_1 + a \rho\_2 t\_2\right) + j \sin\left(a\rho\_1 t\_1 + a \rho\_2 t\_2\right)$$

The 2D NILT function is called from a command line as follows: nilt2c('F',xm,ym); where F is a name of another function in which the *F*(*p,q*) is defined, and xm and ym mark upper limits of the original variables *x* and y. In our case, and for ω1 = ω2 = 2π, this function can have a form

function f=exp2c(p,q) f=1./(p-2\*pi\*j)./(q-2\*pi\*j);

and its calling can look like nilt2c('exp2c',3,3); with graphical results in Fig. 4. As the originals are bounded by values ±1, and α1 = α2 = 0, we can see the errors satisfy (21) very well (δ*M* = 10-8 was considered), excluding beginnings of the 2D region.

0 0.5 1

0 0.5 1

0 0.5 1 1.5

f 4

f 5

f 6

A Matlab language listing is shown in Tab. 3, with all the parameters denoted by similar

Nevertheless, the Laplace transform variables and the original variables have changed in this listing as *s*1 → *p*, *s*2 → *q*, and *t*1 → *x*, *t*2 → *y*, respectively, which simplified a writing. The parameters are then indexed in compliance with these new notations. Besides, numbers of

For the same reasons as at the 1D NILT, the 2D NILT method discussed here enables to treat complex images of two variables resulting in complex originals. We will śhow it on a simple

, = ( , )= *j t t Fs s ft t e sj s j*

( )( ) ( )( )

 ωω

12 1 2 21 12 1 2 2 22 2 2 22 2 1 12 2 1 12 2

− + + ++ ++

ω

( )( )

++ +

The 2D NILT function is called from a command line as follows: nilt2c('F',xm,ym); where F is a name of another function in which the *F*(*p,q*) is defined, and xm and ym mark

and its calling can look like nilt2c('exp2c',3,3); with graphical results in Fig. 4. As

α1 = α

*M* = 10-8 was considered), excluding beginnings of the 2D region.

*ss ss*

 ω 0 10 20 30 -0.5

0 5 10 15 -0.5

0 1 2 3 -0.5

t

0 10 20 30

0 5 10 15

0 1 2 3

t

100 ERRORS

10-10

10-10

10-10

( ) 11 2 2

ω1 = ω2 = 2π

2 = 0, we can see the errors satisfy (21)

. (26)

, this function

+

ω ω

− − . (25)

 ω

ωω

100

100

ORIGINALS

0 1 2 3 4

0 2 4 6 8

0 1 2 3 4

t

( ) ( )( )

After rearranging the above equation, we can also write

( )

, =

ω

( , ) = cos sin

upper limits of the original variables *x* and y. In our case, and for

ωω

Fig. 3. Computed originals and errors for test Laplace transforms in Tab. 2

points used to plot three-dimensional graphical results are set by xpl and ypl.

1 2 1 2 1 12 2 1

1 2 11 2 2 11 2 2

*ft t t t j t t*

*s s s s Fs s j*

ω ω

ωω

100 ERRORS

10-10

10-10

10-10

100

100

0 1 2 3 4 0

0 2 4 6 8 0

0 1 2 3 4 -1

t

way as in the previous case.

ORIGINALS

0.5 1

0.2 0.4

> 0 1

transform pair

can have a form

very well (

function f=exp2c(p,q)

δ

f=1./(p-2\*pi\*j)./(q-2\*pi\*j);

the originals are bounded by values ±1, and

f 3

f 2

f 1

```
% ----- 2D NILT based on partial inversions, by L. Brančík -----
function fxy=nilt2c(F,xm,ym);
alfax=0; alfay=0; Mx=256; My=256; P=3; Er=1e-8; % adjustable
xpl=64; ypl=64; % adjustable
Nx=2*Mx; Ny=2*My; qd=2*P+1; Ke=log(1-1/sqrt(1+Er));
nx=2*xm*Nx/(Nx-2); ny=2*ym*Ny/(Ny-2);
omegax=2*pi/nx; omegay=2*pi/ny; sigx=alfax-Ke/nx;
sigy=alfay-Ke/ny; qd1=qd-1; Nxw=Nx+qd1; Nyw=Ny+qd1;
Asigx=sigx-i*omegax*(0:Nxw); Asigy=sigy-i*omegay*(0:Nyw);
Asigx2=cat(2,Asigx,conj(Asigx));
rx=[1:Mx/xpl:Mx,Mx]; ry=[1:My/ypl:My,My];
x=linspace(0,xm,Mx); y=linspace(0,ym,My); x=x(rx); y=y(ry);
[q,p]=meshgrid(Asigy,Asigx2); Fpq(:,:,1)=feval(F,p,q);
[q,p]=meshgrid(conj(Asigy),Asigx2); Fpq(:,:,2)=feval(F,p,q);
Fpyp=Pnilt(Fpq,Ny,ry,qd,y,ny,omegay,sigy); % Pnilt to get F(p,y)
Fpy(:,:,1)=Fpyp(1:Nxw+1,:).';
Fpy(:,:,2)=Fpyp(Nxw+2:2*Nxw+2,:).';
fxy=Pnilt(Fpy,Nx,rx,qd,x,nx,omegax,sigx); % Pnilt to get f(x,y)
figure; mesh(x,y,real(fxy));
figure; mesh(x,y,imag(fxy)); % optional
% ------ PARTIAL NILT based on FFT/IFFT/Q-D, by L.Brančík ------
function fx=Pnilt(Fq,N,grid,qd,xy,nxy,omega,c);
fx(:,:,1)=fft(Fq(:,:,1),N,2); fx(:,:,2)=N*ifft(Fq(:,:,2),N,2);
fx=fx(:,grid,:); delv=size(Fq,1); delxy=length(xy);
d=zeros(delv,qd,2); e=d; q=Fq(:,N+2:N+qd,:)./Fq(:,N+1:N+qd-1,:);
d(:,1,:)=Fq(:,N+1,:); d(:,2,:)=-q(:,1,:);
for r=2:2:qd-1
 w=qd-r; e(:,1:w,:)=q(:,2:w+1,:)-q(:,1:w,:)+e(:,2:w+1,:);
 d(:,r+1,:)=-e(:,1,:);
 if r>2
 q(:,1:w-1,:)=q(:,2:w,:).*e(:,2:w,:)./e(:,1:w-1,:);
 d(:,r,:)=-q(:,1,:);
 end
end
A2=zeros(delv,delxy,2); B2=ones(delv,delxy,2);
A1=repmat(d(:,1,:),[1,delxy,1]); B1=B2;
z1(1,:,1)=exp(-i*omega*xy); z1(1,:,2)=conj(z1(1,:,1));
z=repmat(z1,[delv,1]);
for n=2:qd
 Dn=repmat(d(:,n,:),[1,delxy,1]);
 A=A1+Dn.*z.*A2; B=B1+Dn.*z.*B2; A2=A1; B2=B1; A1=A; B1=B;
end
fx=fx+A./B; fx=sum(fx,3)-repmat(Fq(:,1),[1,delxy]);
fx=repmat(exp(c*xy)/nxy,[delv,1]).*fx; fx(:,1)=2*fx(:,1);
```
Table 3. Matlab listing of 2D NILT based on partial inversions

Another simple example shows a shifted 2D unit step, with different shifts along the axis. A corresponding transform pair is

$$F(s\_1, s\_2) = \frac{\exp(-2s\_1 - s\_2)}{s\_1 s\_2} \quad \mapsto \quad f(t\_1, t\_2) = \underline{1}(t\_1 - 2, t\_2 - 1) \,. \tag{27}$$

In this case, a displaying imaginary part gives a zero plane, and the respective line in the 2D NILT function can be inactivated. The graphical results are depicted in Fig. 5, including an absolute error. The respective Matlab function can be of a form

Numerical Inverse Laplace Transforms for Electrical Engineering Simulation 61

− − ≈− =