**Dual Heuristic Neural Programming Controller for Synchronous Generator**

Mato Miskovic, Ivan Miskovic and Marija Mirosevic

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/58377

#### **1. Introduction**

24

[22] R. Silva-Ortigoza, C. Márquez-Sánchez, F. Carrizosa-Corral, M. Antonio-Cruz, J. M. Alba-Martínez, and G. Saldaña-González, "Hierarchical velocity control based on differential flatness for a DC/DC Buck converter-DC motor System," *Mathematical Problems in Engineering*, vol. 2014, Article ID 912815, pp. 1–12, 2014. Available at

[23] R. Silva-Ortigoza, V. M. Hernández-Guzmán, M. Antonio-Cruz, and D. Muñoz-Carrillo, "DC/DC Buck power converter as a smooth starter for a DC motor based on a hierarchical control," *IEEE Transactions on Power Electronics*, Article in Press, Available

[24] H. Sira-Ramirez and S. K. Agrawal, *Differentially Flat Systems*, Marcel Dekker, New York,

[26] T. A. Vidal Calleja, *Generalización del método de campos potenciales artificiales para un vehículo articulado*. Tesis de Maestría. Sección de Mecatrónica del Departamento de

Ingeniería Eléctrica del CINVESTAV-IPN, Mexico City, Mexico, 2002.

http://dx.doi.org/10.1155/2014/912815

102 MATLAB Applications for the Practical Engineer

with DOI: 10.1109/TPEL.2014.2311821

http://www.controlautomatico.com.mx/trabajos.html

[25] Supplementary Material. Available at

2004.

This chapter describes the development of voltage control system of a synchronous generator based on neural networks. Recurrent (dynamic) neural networks (RNN) are used, as a type that has great capabilities in approximation of dynamic systems [1]. Two algorithms are used for training – Dual Heuristic Programming (DHP) and Globalized Dual Heuristic Program‐ ming (GDHP). The algorithms have been developed for the optimal control of nonlinear systems using dynamic programming principles.

Neural voltage controller is developed in MATLAB and Simulink environment. For training purposes a mathematical model of synchronous generator is designed and applied in Simu‐ link. DHP and GDHP algorithms are designed in Simulink, with matrix calculations in Sfunctions. Algorithms are used for offline training of neural networks (NN). In the second part, the same functions are redesigned as real time controllers, based on the Real Time Windows Target Toolbox.

DHP or GDHP algorithms provide a significant improvement over conventional PI controller with, in some cases, power system stabilizer (PSS). Conventional linear controller can only provide optimal control in single operating point. Algorithms described in the chapter provide significantly better control over the whole operating range by minimization of the user defined cost function during the training of the neural network.

Digital voltage controller is implemented on real system in two basic steps. First step is neural network design and training in Matlab environment on desktop computer. Second step consists of transfer of controller with computed coefficients to the digital control system hardware.

© 2014 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

For the controller to have optimal performance in all real operating conditions, neural networks must be trained for the whole working range of the plant. Procedure described in the chapter can be applied on standard voltage control systems by replacing the PI controller with the offline-trained neural network. Most modern digital control systems have sufficient hardware and software support for neural network implementation.

#### **2. Recurrent neural networks**

Neural networks are computational models capable of machine learning and pattern recog‐ nition. They can also be used as universal approximators of continuous real functions. Typical dynamic neural network is shown in Figure 1.

Figure 1. Recurrent neural network structure **Figure 1.** Recurrent neural network structure

ϕ is activation function,

activation functions are logsig, tansig and Gaussian.

its derivative () ()2 φ φ u u ′ = −1 is shown in Figure 2.

wkj and x

where

Neural network is formed from interconnected layers of neurons. Output of each neuron is calculated from: Neural network is formed from interconnected layers of neurons. Output of each neuron is calculated from:

kj are weights and outputs from previous layer and b is bias.

Typical activation function maps input values u t() , ∈ −∞ ∞ to output y t( ) 0,1 ∈ or y t( ) 1,1 ∈ − .

Activation function must be derivable over the whole range. Among typically used nonlinear

Dynamic behavior is added by introducing a memory element (step delay) that stores neuron output and reintroduces it as input to all neurons in the same layer in next time step. This type of network is named recurrent or dynamic neural network. An example of tansig

(Hyperbolic Tangent Sigmoid Transfer Function) activation function, ( ) 2 1 2 1

$$\mathbf{y}\_{k} = \phi \begin{pmatrix} \sum \\ j = 0 \end{pmatrix} \begin{pmatrix} \mathbf{w}\_{kj} \mathbf{x}\_{kj} \\ \mathbf{0} \end{pmatrix} + b \begin{pmatrix} \\ \\ \end{pmatrix} \tag{1}$$

u

u

, and

e φ = − − + (1)

where

*wkj*

and *xkj*

derivative *ϕ*(*u*)′

( ) 2

( )

procedures.

networks.

( ) ( )

*u u*

' 1

4 '

*e u e*

1 *u u e*

transfer function

2 1

Hyperbolic tangent sigmoid

Hyperbolic tangent sigmoid transfer derivative function

> - = + = -

*u u*

**Figure 2.** Tansig activation function and its derivative

(1 )

  =1−*ϕ*(*u*)2

*φ* is activation function,

are weights and outputs from previous layer and *b* is bias.

used nonlinear activation functions are *logsig*, *tansig* and Gaussian.

Typical activation function maps input values *u*(*t*)∈ −*∞*, *∞* to output *y*(*t*)∈ 0, 1 or *y*(*t*)∈ −1, 1 . Activation function must be derivable over the whole range. Among typically

Dynamic behavior is added by introducing a memory element (step delay) that stores neuron

This type of network is named recurrent or dynamic neural network. An example of tansig

*u*

Dual Heuristic Neural Programming Controller for Synchronous Generator

*u*

1 + *e* −2*u* −1, and its

http://dx.doi.org/10.5772/58377

105

output and reintroduces it as input to all neurons in the same layer in next time step.

(*Hyperbolic Tangent Sigmoid Transfer Function*) activation function, *ϕ*(*u*)= 2

is shown in Figure 2.

(*u*)

'(*u*) 0.2 0.4 0.6 0.8 1

Input values are in range of *u* ∈ −*∞*, *∞* , with output values in range *φ*(*u*)∈ −1, 1 .

Use of function with simple derivative significantly improves calculation times of numerical

Dynamic (recurrent) neural networks provide good approximation of time-dependent functions. However, training of recurrent neural networks is more complex than training static where

For the controller to have optimal performance in all real operating conditions, neural networks must be trained for the whole working range of the plant. Procedure described in the chapter can be applied on standard voltage control systems by replacing the PI controller with the offline-trained neural network. Most modern digital control systems have sufficient

Neural networks are computational models capable of machine learning and pattern recog‐ nition. They can also be used as universal approximators of continuous real functions. Typical

> ij wh

ΣΣ ϕ

ΣΣ ϕ

ΣΣ ϕ

Figure 1. Recurrent neural network structure

( 1) ij xt w −

Neural network is formed from interconnected layers of neurons. Output of each neuron is

Neural network is formed from interconnected layers of neurons. Output of each neuron is

0

y wx b k kj kj j

æ ö = + ç ÷ å è ø =

( ) 0

Activation function must be derivable over the whole range. Among typically used nonlinear

Dynamic behavior is added by introducing a memory element (step delay) that stores neuron output and reintroduces it as input to all neurons in the same layer in next time step. This type of network is named recurrent or dynamic neural network. An example of tansig

(Hyperbolic Tangent Sigmoid Transfer Function) activation function, ( ) 2 1 2 1

= + ∑

=

kj are weights and outputs from previous layer and b is bias.

to output y t( ) 0,1 ∈

1 Y t( ) x t( )

ij

wo

ΣΣ ϕ

u

(1)

m

*m y wx b k kj kj j*

ϕ

j

hardware and software support for neural network implementation.

**2. Recurrent neural networks**

104 MATLAB Applications for the Practical Engineer

1 u t( )

2 u t( )

calculated from:

ϕ is activation function,

calculated from:

Typical activation function maps input values u t() , ∈ −∞ ∞

**Figure 1.** Recurrent neural network structure

activation functions are logsig, tansig and Gaussian.

its derivative () ()2 φ φ u u ′ = −1 is shown in Figure 2.

wkj and x

where

dynamic neural network is shown in Figure 1.

ij wh

*φ* is activation function,

*wkj* and *xkj* are weights and outputs from previous layer and *b* is bias.

Typical activation function maps input values *u*(*t*)∈ −*∞*, *∞* to output *y*(*t*)∈ 0, 1 or *y*(*t*)∈ −1, 1 . Activation function must be derivable over the whole range. Among typically used nonlinear activation functions are *logsig*, *tansig* and Gaussian.

Dynamic behavior is added by introducing a memory element (step delay) that stores neuron output and reintroduces it as input to all neurons in the same layer in next time step.

This type of network is named recurrent or dynamic neural network. An example of tansig (*Hyperbolic Tangent Sigmoid Transfer Function*) activation function, *ϕ*(*u*)= 2 1 + *e* −2*u* −1, and its derivative *ϕ*(*u*)′ =1−*ϕ*(*u*)2 is shown in Figure 2.

#### **Figure 2.** Tansig activation function and its derivative

(1)

or y t( ) 1,1 ∈ − .

u

, and

e φ = − − + Input values are in range of *u* ∈ −*∞*, *∞* , with output values in range *φ*(*u*)∈ −1, 1 .

Use of function with simple derivative significantly improves calculation times of numerical procedures.

Dynamic (recurrent) neural networks provide good approximation of time-dependent functions. However, training of recurrent neural networks is more complex than training static networks.

#### **2.1. Recurrent neural networks training**

There are many methods for training neural networks, and most of them rely on some form of backpropagation – calculation of partial derivative of system output over network weight coefficients ( ∂ *y* ∂*w* ). Partial derivatives are determined using chain rule, by finding derivatives of weight coefficients in output layer, using intermediate results to calculate next layer, with respect to network structure. Calculation is repeated until all layers are processed.

Described procedure is valid for static neural networks, but is of no use in recurrent networks, as each neuron has inputs from previous time step from the same layer. Two procedures are often used in calculation of partial derivatives: *Backpropagation Through Time* (*BPTT*) and *Real Time Recurrent Learning* (*RTRL*).

In BPTT recurrent neuron inputs are replaced with the same neural network delayed by one time step. Process is iterated N times. Thus, recurrent network is approximated by N feedfor‐ ward networks with delayed outputs. Network weight coefficients, recurrent inputs and outputs must be preserved for past N time steps. Obtained derivative ( ∂ *y*(*t*) ∂*w*(*t*) ) is used to determine new weight coefficients. One advantage of BPTT is that all methods used for training feedforward networks can be utilized in training recurrent networks.

This chapter will describe use of RTRL procedure in training dynamic networks. Output of recurrent network with one hidden and one output layer has the following form:

$$\begin{aligned} \mathbf{x}(t) &= \tanh\left(W\_{\mathbf{l}} \cdot p\_{\mathbf{l}}\right) \\ \mathbf{y}(t) &= W\_{\mathbf{2}} \cdot \mathbf{x}(t) \end{aligned} \tag{2}$$

For hidden layer, partial derivatives are

*NW* 1– number of neurons in hidden layer,

*Nin*– number of neurons in input layer.

**2.2. Neural network training algorithms**

between real and desired network output can be expressed as

Differentiation of (5) provides update values of weight coefficients:

Parameter *β* defines learning rate and must be in range 0<*β* <1.

where

Value of

∂ *xl* (*t*) ∂*wi*, *j*

*2.2.1. Gradient descent*

*NW* 1, *NW* 1 ⋅ (*Nin* + *NW* 1 + 1) .

2 2 () () () () () () ¶ ¶¶ = ¶ ¶¶ *h yt yt xt*

,

*in*

*w t w t* (5)

 + =

Partial derivative of hidden layer output over weight coefficients is

[ ] 1

, , 1 ( ) ( 1) ()' () () ( ) ( 1) d

is calculated using four nested loops (Appendix A).

Matrix of output derivatives over weight coefficients is of the form

Gradient descent is the most commonly used method in neural networks training. Difference

ˆ( ) ( ) ( 1) ( ) b

*ij*

¶ = -- ¶ *ij ij*

*y t wt wt*

*l k i j li lN k i j k i j x t x t f xt p t wt*

¶ æ ö ¶ - = + ç ÷ ¶ ¶ - è ø å *W*

*N*

*wt wt wt* (4)

Dual Heuristic Neural Programming Controller for Synchronous Generator

http://dx.doi.org/10.5772/58377

107

*et yt yt* () () () = - ˆ (6)

*w t* (7)

where

*p*1– hidden layer input, *p*1 = *u*(*t*) ; *x*(*t* −1) ; 1 ,

*x*(*t*)– hidden layer output,

y(t) – network output,

*W*1– hidden layer weight coefficients,

*W*2– output layer weight coefficients, *p*2 = *x*(*t*) ; 1

Equation (2) can easily be expanded for multi layer recurrent networks. Partial derivative of network output over output layer weight coefficients is

$$\frac{\partial \hat{\boldsymbol{y}}(t)}{\partial \boldsymbol{w}\_{y}^{o}(t)} = p\_{2}(t) \tag{3}$$

For hidden layer, partial derivatives are

$$\frac{\partial \mathbf{\hat{y}}(t)}{\partial \mathbf{w}\_h(t)} = \frac{\partial \mathbf{y}(t)}{\partial w\_2(t)} \frac{\partial \mathbf{x}(t)}{\partial w\_2(t)}\tag{4}$$

Partial derivative of hidden layer output over weight coefficients is

$$\frac{\partial \mathbf{\hat{x}}\_{l}(t)}{\partial \mathbf{w}\_{l,j}(t)} = f\left[\mathbf{x}(t)\right]!\_{l} \left(p\_{j}(t)\delta\_{ll} + \sum\_{k=1}^{N\_{W1}} \mathbf{w}(t)\_{l,N\_{k}+k} \frac{\partial \mathbf{x}\_{k}(t-l)}{\partial w\_{l,j}(t-l)}\right) \tag{5}$$

where

**2.1. Recurrent neural networks training**

106 MATLAB Applications for the Practical Engineer

*Time Recurrent Learning* (*RTRL*).

*p*1– hidden layer input, *p*1 = *u*(*t*) ; *x*(*t* −1) ; 1 ,

*W*2– output layer weight coefficients, *p*2 = *x*(*t*) ; 1

network output over output layer weight coefficients is

*W*1– hidden layer weight coefficients,

*x*(*t*)– hidden layer output,

y(t) – network output,

coefficients ( ∂ *y*

where

There are many methods for training neural networks, and most of them rely on some form of backpropagation – calculation of partial derivative of system output over network weight

of weight coefficients in output layer, using intermediate results to calculate next layer, with

Described procedure is valid for static neural networks, but is of no use in recurrent networks, as each neuron has inputs from previous time step from the same layer. Two procedures are often used in calculation of partial derivatives: *Backpropagation Through Time* (*BPTT*) and *Real*

In BPTT recurrent neuron inputs are replaced with the same neural network delayed by one time step. Process is iterated N times. Thus, recurrent network is approximated by N feedfor‐ ward networks with delayed outputs. Network weight coefficients, recurrent inputs and

determine new weight coefficients. One advantage of BPTT is that all methods used for training

This chapter will describe use of RTRL procedure in training dynamic networks. Output of

Equation (2) can easily be expanded for multi layer recurrent networks. Partial derivative of

2 ( ) ( ) ( ) ¶ = ¶ *o ij*

*y t p t*

*yt W xt* (2)

*w t* (3)

respect to network structure. Calculation is repeated until all layers are processed.

outputs must be preserved for past N time steps. Obtained derivative ( ∂ *y*(*t*)

recurrent network with one hidden and one output layer has the following form:

( ) tanh( ) 1 1 () () 2 = × = × *xt W p*

feedforward networks can be utilized in training recurrent networks.

∂*w* ). Partial derivatives are determined using chain rule, by finding derivatives

∂*w*(*t*) ) is used to

*NW* 1– number of neurons in hidden layer,

*Nin*– number of neurons in input layer.

Value of ∂ *xl* (*t*) ∂*wi*, *j* is calculated using four nested loops (Appendix A).

Matrix of output derivatives over weight coefficients is of the form *NW* 1, *NW* 1 ⋅ (*Nin* + *NW* 1 + 1) .

#### **2.2. Neural network training algorithms**

#### *2.2.1. Gradient descent*

Gradient descent is the most commonly used method in neural networks training. Difference between real and desired network output can be expressed as

$$e(t) = \mathbf{y}(t) - \hat{\mathbf{y}}(t) \tag{6}$$

Differentiation of (5) provides update values of weight coefficients:

$$\mathcal{W}\_{\mathcal{Y}}(t) = \mathcal{W}\_{\mathcal{Y}}(t-\mathrm{l}) - \beta \frac{\widehat{\mathcal{O}}\widehat{\mathcal{Y}}(t)}{\widehat{\mathcal{O}}\mathcal{W}\_{\mathcal{Y}}(t)} \tag{7}$$

Parameter *β* defines learning rate and must be in range 0<*β* <1.

#### *2.2.2. Kalman filter based training*

Kalman filter [5] was developed as a recursive procedure of optimal filtering of linear dynamic state-space systems. For RNN training extended Kalman filter (EKF) is used. It is also appli‐ cable to nonlinear systems using linearization around operating point. The algorithm is used as a predictive-corrective procedure.

Kalman filter based recurrent network training is described in [6] and [7].

Training is performed in two steps. The first step is calculation of prediction values of weight coefficients. Second step is correction of predicted variables based on measured state space.

Algorithm is defined by the following relations:

$$K(t) = P(t)H(t)\left(\eta(t)I + H^T(t)P(t)H(t)\right) \tag{8}$$

$$W(t) = W(t-\mathbf{l}) + K(t)\varepsilon(t) \tag{9}$$

[ ] [ ] 1 () () () () () ˆ ˆ 2 =× - × - *T Et yt yt yt yt* (11)

Dual Heuristic Neural Programming Controller for Synchronous Generator

system

*u t*( ) *y t*( )

http://dx.doi.org/10.5772/58377

109

1 *z*-

*y t* ˆ( )

Neural Networks

1 *z*-

Er

ror

^(*t*) is the difference of the measured system outputs and output values of the

Dynamic system identification can be performed according to [9], using either series-parallel

Series- parallel Identification Parallel Identification

In series-parallel structure Figure 3, model neural network input vector is expanded with system outputs from previous steps. In parallel structure Figure 3, network output from previous states is added to input vector. Series-parallel identification provides estimated system output for next time step only, while parallel identification can predict multiple future iterations. However, parallel structure is harder to train and is susceptible to state drifting.

To verify the validity of the RTRL algorithm, neural networks were trained on several discrete mathematical functions. Supervised training has been used, with control value expressed as

A parallel identification structure was used on dynamic SISO system

four neurons in hidden layer, and one neuron in linear output layer. Identification results are

. Training has been performed on Recurrent Neural Network with

Selected cost function represents squared error of the identified model.

*u t*( ) *y t*( ) Non linear dynamic

1 *z*-

*y t* ˆ( )

where *y*(*t*)− *y*

or parallel structure.

1 *z*-

*y*(*t* + 1)= *y*(*t*)

shown in Figure 4.

1 + *y*(*t*)

2 + *u*(*t*) 3

Non linear dynamic system

Neural Networks

**Figure 3.** Nonlinear dynamic system identification

Error

difference between desired and actual network output.

model.

$$P(t+l) = P(t) - K(t)H^T(t)P(t) + Q(t) \tag{10}$$

where

*K*(*t*)– Kalman gain matrix, determined form network output derivatives over weight coeffi‐ cients *H* (*t*)= ∂ *y*(*t*) ∂*w* . System linearization is accomplished by calculating *H* (*t*).

*P*(*t*)– Error covariance matrix, calculated in every iteration

To escape local minima during training, a diagonal matrix *Q*(*t*) is added to covariance matrix *P*(*t*) in each time step.

Updates to weight coefficients are determined from (9), as a product of Kalman gain and *ε*(*t*) in current step.

Training based on Kalman filter is fast, has good rate of convergence and is often the only practical solution for recurrent networks training. However, the procedure is demanding in terms of computational power and memory requirements.

#### **3. System identification**

Nonlinear dynamic systems are identified as input-output system model, or as a state space system. For the selected process, identification is performed by minimization of the cost function

$$E(t) = \frac{1}{2} \cdot \left[\boldsymbol{\chi}(t) - \boldsymbol{\hat{\chi}}(t)\right]^T \cdot \left[\boldsymbol{\chi}(t) - \boldsymbol{\hat{\chi}}(t)\right] \tag{11}$$

where *y*(*t*)− *y* ^(*t*) is the difference of the measured system outputs and output values of the model.

Selected cost function represents squared error of the identified model.

Dynamic system identification can be performed according to [9], using either series-parallel or parallel structure.

*2.2.2. Kalman filter based training*

108 MATLAB Applications for the Practical Engineer

as a predictive-corrective procedure.

where

cients *H* (*t*)= ∂ *y*(*t*)

*P*(*t*) in each time step.

**3. System identification**

in current step.

function

Algorithm is defined by the following relations:

*P*(*t*)– Error covariance matrix, calculated in every iteration

terms of computational power and memory requirements.

Kalman filter [5] was developed as a recursive procedure of optimal filtering of linear dynamic state-space systems. For RNN training extended Kalman filter (EKF) is used. It is also appli‐ cable to nonlinear systems using linearization around operating point. The algorithm is used

Training is performed in two steps. The first step is calculation of prediction values of weight coefficients. Second step is correction of predicted variables based on measured state space.

> () () () () () () () = + (h

> > *Wt Wt Kt t* ( ) ( 1) ( ) ( ) = -+

*K*(*t*)– Kalman gain matrix, determined form network output derivatives over weight coeffi‐

∂*w* . System linearization is accomplished by calculating *H* (*t*).

To escape local minima during training, a diagonal matrix *Q*(*t*) is added to covariance matrix

Updates to weight coefficients are determined from (9), as a product of Kalman gain and *ε*(*t*)

Training based on Kalman filter is fast, has good rate of convergence and is often the only practical solution for recurrent networks training. However, the procedure is demanding in

Nonlinear dynamic systems are identified as input-output system model, or as a state space system. For the selected process, identification is performed by minimization of the cost

e

) *T Kt PtHt tI H tPtHt* (8)

( 1) ( ) ( ) ( ) ( ) ( ) += - + *T Pt Pt KtH tPt Qt* (10)

(9)

Kalman filter based recurrent network training is described in [6] and [7].

In series-parallel structure Figure 3, model neural network input vector is expanded with system outputs from previous steps. In parallel structure Figure 3, network output from previous states is added to input vector. Series-parallel identification provides estimated system output for next time step only, while parallel identification can predict multiple future iterations. However, parallel structure is harder to train and is susceptible to state drifting.

To verify the validity of the RTRL algorithm, neural networks were trained on several discrete mathematical functions. Supervised training has been used, with control value expressed as difference between desired and actual network output.

A parallel identification structure was used on dynamic SISO system *y*(*t* + 1)= *y*(*t*) 1 + *y*(*t*) 2 + *u*(*t*) 3 . Training has been performed on Recurrent Neural Network with four neurons in hidden layer, and one neuron in linear output layer. Identification results are shown in Figure 4.

0 50 100 150 200 250

0 50 100 150 200 250

0 0.5 1 1.5 2

(c)

http://dx.doi.org/10.5772/58377

111

k

(*t*) ). After finding output derivatives, various

*u t* (12)

Dual Heuristic Neural Programming Controller for Synchronous Generator

E(t)

u(t) y(t) y(t) RNN

x 104

u(t) y(t) y(t) RNN

k

Both examples were designed in Matlab and Simulink. Neural network was trained using clanguage S-function (Appendix B). Although the training took 10000 steps, all important

One of the most important features of neural network modeling is the possibility to determine

(*t*) ∂*uj*

( ) 1 , ,

1 1 ( ) '( ) ( ) = = æ ö ¶ =ç × × ÷ ¶ ç ÷ è ø å å *NW Np i ij k kj*

*y t W f xW*

*j j k*

Input signal used for training and verification is the same as in previous example:

22 *kTD*), where (*TD* =1*s*).

**3.1. Determination of partial derivatives of functions using neural networks**

system dynamics were visible on trained network after 2000 steps.

(b)

k

(a)

2*π*

**Figure 5.** Identification of dynamic function 2

50 *kTD*) + sin(

2*π*

system output derivatives over inputs ( ∂ *yi*

optimization algorithms can be used.

Partial derivative is calculated from (2)

u(t)

u(t)

*u*(*t*)=sin(

k(t)

y(t)

**Figure 4.** Identification of dynamic function 1

Figure 4a compares initial response of model output *y*(*t*) and network output *y* ^(*t*). Figure 4b represents the same values with trained network. There is a perfect match between function output values and output of the trained NN. Figure 4c displays training rate as error square (*y*(*t*)− *y* ^(*t*))2 . Training was completed in 20000 steps.

Input signal used for training and verification is *u*(*t*)=sin( 2*π* 50 ) + sin( 2*π* 12 ), with sample time (*TD* =1*s*).

In second test, a system with one input, two state variables and one output has been identified:

$$\begin{aligned} \mathbf{x}\_1(t+1) &= 0.5\mathbf{x}\_2 \cdot (t) + 0.2 \cdot \mathbf{x}\_1(t) \cdot \mathbf{x}\_2(t) \\ \mathbf{x}\_2(t+1) &= -0.3 \cdot \mathbf{x}\_1(t) + 0.8 \cdot \mathbf{x}\_2(t) + \boldsymbol{\mu}(t) \\ \mathbf{y}(t+1) &= (\mathbf{x}\_1(t) + \mathbf{x}\_2(t))^2 \end{aligned}$$

System was identified using series-parallel structure with RNN with 5 neurons in hidden layer, and one linear neuron in output layer. Identification results are shown in Figure 5.

**Figure 5.** Identification of dynamic function 2

0 20 40 60 80 100 120 140 160 180 200

0 20 40 60 80 100 120 140 160 180 200 -100

represents the same values with trained network. There is a perfect match between function output values and output of the trained NN. Figure 4c displays training rate as error square

In second test, a system with one input, two state variables and one output has been identified:

System was identified using series-parallel structure with RNN with 5 neurons in hidden layer,

and one linear neuron in output layer. Identification results are shown in Figure 5.

Figure 4a compares initial response of model output *y*(*t*) and network output *y*

. Training was completed in 20000 steps.

Input signal used for training and verification is *u*(*t*)=sin(

(b)

k

u(t) y(u) y(u) RNN

0 0.5 1 1.5 2

(c)

k

2*π* 50 ) + sin( 2*π*

E(t)

x 104

^(*t*). Figure 4b

12 ), with sample time

u(t) y(t) y(t) RNN

k

(a)

110 MATLAB Applications for the Practical Engineer

**Figure 4.** Identification of dynamic function 1

*x*1(*t* + 1)=0.5*x*2 ⋅ (*t*) + 0.2⋅ *x*1(*t*) ⋅ *x*2(*t*) *x*2(*t* + 1)= −0.3⋅ *x*1(*t*) + 0.8⋅ *x*2(*t*) + *u*(*t*)

0

u(t)

(*y*(*t*)− *y*

(*TD* =1*s*).

^(*t*))2

*y*(*t* + 1)=(*x*1(*t*) + *x*2(*t*))2

y(t)

50

100

u(t)

y(t)

Input signal used for training and verification is the same as in previous example:

$$u(t) = \sin(\frac{2\pi}{50}kT\_D) + \sin(\frac{2\pi}{22}kT\_D), \text{ where } \{T\_D = 1s\}.$$

Both examples were designed in Matlab and Simulink. Neural network was trained using clanguage S-function (Appendix B). Although the training took 10000 steps, all important system dynamics were visible on trained network after 2000 steps.

#### **3.1. Determination of partial derivatives of functions using neural networks**

One of the most important features of neural network modeling is the possibility to determine system output derivatives over inputs ( ∂ *yi* (*t*) ∂*uj* (*t*) ). After finding output derivatives, various optimization algorithms can be used.

Partial derivative is calculated from (2)

$$\frac{\partial \mathbf{y}\_t(t)}{\partial u\_f(t)} = \sum\_{f=1}^{N\_{W1}} \left( W\_{l,f} \cdot \sum\_{k=1}^{N\_p} \left( f\_k \, \mathbf{y}(\mathbf{x}) \cdot W\_{k,f} \right) \right) \tag{12}$$

where *Np* = *NW* 1 + *Nin* + 1 is number of hidden layer inputs. Partial derivative ( ∂ *yi* (*t*) ∂*uj* (*t*) ) is obtained from network backpropagation.

neural networks and using them in the structure according to the selected ACD algorithm. ACD algorithms are implemented using a heuristic approach – system state space variables, control values and cost function are selected based on process characteristics. This text focuses on the Heuristic Dynamic Programming (HDP), Dual Heuristic Programming (DHP) and

Heuristic dynamic programming (*HDP*) is the simplest (by structure) ACD algorithm. The basis of all HDP algorithms is minimization of Bellman dynamic programming function *J*(*t*):

*γ* represents discount factor with value in range 0<*γ* <1, ensuring the convergence of Bellman

Control signal that minimizes Bellman sum is generated using one of the Adaptive Critic algorithms. In practical control systems, HDP optimization is realized with three neural

Model neural network in *7a* is used for model identification. Network inputs are control signal *u*(*t*) and output state in the previous time step *y*(*t* −1). Network output represents estimated

to be estimated. Model NN is structured according to the selected process model, where the control value is input to the NN, and NN output represents predicted system output. The same

> = × å () () *M T M M t*

For adaptive control systems, model NN is trained continuously, allowing for adaptation to change of system parameters, and optimal control. Model neural network is trained simulta‐

*Jt Ut k* (13)

Dual Heuristic Neural Programming Controller for Synchronous Generator

http://dx.doi.org/10.5772/58377

113

^(*t*). Only systems outputs used in utility function need

*E E tE t* (14)

= - () () ˆ *E yt yt M* (15)

0 () ( ) g

*U* (*t* + 1) is user defined utility function based on system state space variables.

¥ = = + å *k k*

General Dual Heuristic Programming (GDHP) [1,2].

**4.1. Heuristic dynamic programming**

sum *J*(*t*) and cancellation of old values,

model output in the current time step *y*

neously with other networks over time.

networks with a structure as shown in Figure 7.

NN output values are used to calculate the cost function.

Model neural network is trained to minimize the quadratic criterion:

Critic neural network in *7b* is used for identification of Bellman sum *J*(*t*).

where

where:

For nonlinear function *y*(*t*)=*u*(*t* −1) 2 + 6⋅*u*(*t* −1) 3 , partial derivative ∂ *y*(*t*) ∂*u*(*t*) was first calculated analytically, and then obtained from model network using (12). Results are compared in Figure 6a and 6b.

**Figure 6.** Identification of nonlinear function and calculation of partial derivatives using neural networks

Results from Figure 6b show good match between calculated function derivatives and values obtained from neural network of identified model.

In S-function for training RNN (Appendix B) a procedure is shown for online calculation of partial derivatives of ( ∂ *yi* (*t*) ∂*uj* (*t*) ), using (12).

#### **4. Adaptive critic design**

Adaptive critic design (ACD) is a set of optimal control procedures of nonlinear systems. It is based on dynamic programming and neural networks. Optimal control is achieved by training neural networks and using them in the structure according to the selected ACD algorithm. ACD algorithms are implemented using a heuristic approach – system state space variables, control values and cost function are selected based on process characteristics. This text focuses on the Heuristic Dynamic Programming (HDP), Dual Heuristic Programming (DHP) and General Dual Heuristic Programming (GDHP) [1,2].

#### **4.1. Heuristic dynamic programming**

Heuristic dynamic programming (*HDP*) is the simplest (by structure) ACD algorithm. The basis of all HDP algorithms is minimization of Bellman dynamic programming function *J*(*t*):

$$J(t) = \sum\_{k=0}^{n} \gamma^k U(t+k) \tag{13}$$

where

where *Np* = *NW* 1 + *Nin* + 1 is number of hidden layer inputs. Partial derivative ( ∂ *yi*

2 + 6⋅*u*(*t* −1)

3

3 3.01 3.02 3.03 3.04 3.05 3.06 3.07 3.08 3.09 3.1

3 3.01 3.02 3.03 3.04 3.05 3.06 3.07 3.08 3.09 3.1

k

Results from Figure 6b show good match between calculated function derivatives and values

In S-function for training RNN (Appendix B) a procedure is shown for online calculation of

Adaptive critic design (ACD) is a set of optimal control procedures of nonlinear systems. It is based on dynamic programming and neural networks. Optimal control is achieved by training

**Figure 6.** Identification of nonlinear function and calculation of partial derivatives using neural networks

obtained from neural network of identified model.

(*t*) ∂*uj*

(*t*) ), using (12).

k

(b)

analytically, and then obtained from model network using (12). Results are compared in Figure

(a)

, partial derivative ∂ *y*(*t*)

u(k) Y(k) y(k) NN

obtained from network backpropagation.

For nonlinear function *y*(*t*)=*u*(*t* −1)

112 MATLAB Applications for the Practical Engineer

6a and 6b.

partial derivatives of ( ∂ *yi*

**4. Adaptive critic design**

dY/du

u(k)

Y(k)

(*t*) ∂*uj*

∂*u*(*t*) was first calculated

x 104

dY/du dY/du RNN

x 104

(*t*) ) is

*γ* represents discount factor with value in range 0<*γ* <1, ensuring the convergence of Bellman sum *J*(*t*) and cancellation of old values,

*U* (*t* + 1) is user defined utility function based on system state space variables.

Control signal that minimizes Bellman sum is generated using one of the Adaptive Critic algorithms. In practical control systems, HDP optimization is realized with three neural networks with a structure as shown in Figure 7.

Model neural network in *7a* is used for model identification. Network inputs are control signal *u*(*t*) and output state in the previous time step *y*(*t* −1). Network output represents estimated model output in the current time step *y* ^(*t*). Only systems outputs used in utility function need to be estimated. Model NN is structured according to the selected process model, where the control value is input to the NN, and NN output represents predicted system output. The same NN output values are used to calculate the cost function.

Model neural network is trained to minimize the quadratic criterion:

$$\left\| E\_M \right\| = \sum\_t E\_M^T(t) \cdot E\_M(t) \tag{14}$$

where:

$$E\_M = \mathbf{y}(t) - \hat{\mathbf{y}}(t) \tag{15}$$

For adaptive control systems, model NN is trained continuously, allowing for adaptation to change of system parameters, and optimal control. Model neural network is trained simulta‐ neously with other networks over time.

Critic neural network in *7b* is used for identification of Bellman sum *J*(*t*).

where

4.1 Heuristic Dynamic Programming

(a) (b)

Heuristic dynamic programming (HDP) is the simplest (by structure) ACD algorithm. The basis of all HDP algorithms is minimization of Bellman dynamic programming function J t( ) :

> 0 () ( ) k k Jt Ut k γ ∞ =

γ represents discount factor with value in range 0 1 < < γ , ensuring the convergence of Bellman sum J t( ) and cancellation of old values, U t( 1) + is user defined utility function based on system state space variables. Control signal that minimizes Bellman sum is generated using one of the Adaptive Critic

= + ∑ (13)

where:

Partial derivatives ∂ *y*(*t*)

**4.2. Dual heuristic programming**

of critic network.

with

where ∂ *J y*

^(*t*) ∂ *y*

procedure efficiency by increasing training speed.

Structure of DHP control system is shown in Figure 8.

Training of critic neural network minimizes scalar product:

∂*u*(*t*) and ∂ *J*(*t*)

gation in the manner described in section 2.5.

() () ( ) () () ¶ ¶ = × ¶ ¶ *A yt Jt E t*

Minimization of scalar product *EA* in (18) by training of action network results in minimi‐

Action neural network is trained from the output network (unsupervised algorithms). Training of neural networks is carried out using (19) with the training concept shown in Figure 7c.

Dual Heuristic Programming optimal control algorithm is developed with a goal of increasing

Model network training is the same as in HDP procedure. The improvement is in the direct

= × å () () *T C CC t*

[ ] ˆ ˆ ( ) ( 1) ( ) [ ] [ ] ( ) ˆ() () () g

^(*t*) represents partial derivative of Bellman sum over state space variables, and

¶ ¶ +¶ =- - ¶ ¶¶ *C J yt J yt U yt E t*

is estimated from critic network. The second part of (21) is determined by deriving

calculation of the Bellman sum partial derivative over state space variables ∂ *J*(*t* + 1)

zation of Bellman sum in *J*(*t*) in (13). Minimum of *J*(*t*) is achieved when ∂ *J*(*t*)

network training criterion approaches zero, thus ensuring end of training process.

*ut yt* (19)

Dual Heuristic Neural Programming Controller for Synchronous Generator

∂ *y*(*t*) ≈0, as action

http://dx.doi.org/10.5772/58377

115

∂ *y*(*t*) as output

∂ *y*(*t*) are determined from model network using backpropa‐

*E EtEt* (20)

*yt yt yt* (21)

Figure 7. HDP Control structure

(c)

model output in the current time step y t ˆ( ) . Only systems outputs used in utility function need

**Figure 7.** HDP Control structure

Training is used to minimize Model neural network in 7a is used for model identification. Network inputs are control signal u t( ) and output state in the previous time step y t( 1) − . Network output represents estimated

$$\left\|E\_C\right\| = \sum\_t E\_C^T(t) \cdot E\_C(t) \tag{16}$$

where:

$$E\_C(t) = J(t) - \left[\mathcal{Y} \cdot J(t+1) + U(t)\right] \tag{17}$$

Critic neural network training is performed with a copy of the same network which is fed with predicted future system outputs *y* ^(*t* + 1), thus approximating *J*(*t* + 1).

Action neural network training is based on model and critic network outputs, by minimizing criterion:

$$\left\| E\_{\mathcal{A}} \right\| = \sum\_{t} E\_{\mathcal{A}}^{T}(t) \cdot E\_{\mathcal{A}}(t) \tag{18}$$

where:

$$E\_{\mathcal{A}}(t) = \frac{\hat{\mathcal{O}}\mathcal{Y}(t)}{\hat{\mathcal{O}}u(t)} \cdot \frac{\hat{\mathcal{O}}J(t)}{\hat{\mathcal{O}}\mathcal{Y}(t)}\tag{19}$$

Partial derivatives ∂ *y*(*t*) ∂*u*(*t*) and ∂ *J*(*t*) ∂ *y*(*t*) are determined from model network using backpropa‐ gation in the manner described in section 2.5.

Minimization of scalar product *EA* in (18) by training of action network results in minimi‐ zation of Bellman sum in *J*(*t*) in (13). Minimum of *J*(*t*) is achieved when ∂ *J*(*t*) ∂ *y*(*t*) ≈0, as action network training criterion approaches zero, thus ensuring end of training process.

Action neural network is trained from the output network (unsupervised algorithms). Training of neural networks is carried out using (19) with the training concept shown in Figure 7c.

#### **4.2. Dual heuristic programming**

Dual Heuristic Programming optimal control algorithm is developed with a goal of increasing procedure efficiency by increasing training speed.

Model network training is the same as in HDP procedure. The improvement is in the direct calculation of the Bellman sum partial derivative over state space variables ∂ *J*(*t* + 1) ∂ *y*(*t*) as output of critic network.

Structure of DHP control system is shown in Figure 8.

Training of critic neural network minimizes scalar product:

$$\left\| E\_C \right\| = \sum\_t E\_C^T(t) \cdot E\_C(t) \tag{20}$$

with

Training is used to minimize

**Figure 7.** HDP Control structure

4.1 Heuristic Dynamic Programming

networks with a structure as shown in Figure 7.

MODEL NETWORK

y(t)

NETWORK

CRITIC

where

y t( 1) −

TD

TD

114 MATLAB Applications for the Practical Engineer

TD

TD

TD

TD

TD

y t( 1) −

TD

TD

TD

u t( )

u t( )

Heuristic dynamic programming (HDP) is the simplest (by structure) ACD algorithm. The basis of all HDP algorithms is minimization of Bellman dynamic programming function J t( ) :

> 0 () ( ) k k Jt Ut k γ ∞ =

γ represents discount factor with value in range 0 1 < < γ , ensuring the convergence of Bellman sum J t( ) and cancellation of old values, U t( 1) + is user defined utility function based on system state space variables. Control signal that minimizes Bellman sum is generated using one of the Adaptive Critic algorithms. In practical control systems, HDP optimization is realized with three neural

y t ˆ( )

yt yt () () − ˆ y t( )

= + ∑ (13)

Jt Jt Ut ( ) ( 1) ( ) − ⋅ ++ [γ ]

CRITIC

ACTION

CRITIC

J t( 1) +

J t( )

predicted future system outputs *y*

where:

criterion:

= × å () () *T C CC t*

Figure 7. HDP Control structure

J t( )

(c)

( ) ( ) y t u t ∂ ∂

( ) ( ) J t y t ∂ ∂

u t ˆ( ) MODEL y t ˆ( )

TD

TD

(a) (b)

y t( )

X

TD

y t ˆ( 1) +

TD

TD

TD

y t ˆ( )

Model neural network in 7a is used for model identification. Network inputs are control signal u t( ) and output state in the previous time step y t( 1) − . Network output represents estimated model output in the current time step y t ˆ( ) . Only systems outputs used in utility function need

> *E t Jt Jt Ut C* ( ) ( ) ( 1) ( ) = - × ++ [g

Critic neural network training is performed with a copy of the same network which is fed with

Action neural network training is based on model and critic network outputs, by minimizing

= × å () () *T A AA t*

^(*t* + 1), thus approximating *J*(*t* + 1).

*E EtEt* (16)

*E EtEt* (18)

] (17)

$$E\_C(t) = \frac{\hat{\varepsilon}J\left[\hat{\jmath}(t)\right]}{\hat{\varepsilon}\hat{\jmath}(t)} - \gamma \frac{\hat{\varepsilon}J\left[\hat{\jmath}(t+l)\right]}{\hat{\varepsilon}\wp(t)} - \frac{\hat{\varepsilon}U\left[\wp(t)\right]}{\hat{\varepsilon}\wp(t)}\tag{21}$$

where ∂ *J y* ^(*t*) ∂ *y* ^(*t*) represents partial derivative of Bellman sum over state space variables, and is estimated from critic network. The second part of (21) is determined by deriving

4.2 Dual Heuristic Programming

procedure efficiency by increasing training speed.

Structure of DHP control system is shown in Figure 8.

(a)

Dual Heuristic Programming optimal control algorithm is developed with a goal of increasing

( ) J t y t ∂ + ∂

as

Third part of equation (22) is determined from deriving

from action network. Prediction of states *y*

Thus, equation (22) can be transformed to

g

Training of Action Network minimizes scalar product:

∂ *y* ^ *i* (*t* + 1)

∂*u* ^

network.

∂ *J y* ^(*t*) ∂ *y* ^(*t*) ≈0.

with

linear controller.

**4.3. Global dual heuristic programming**

*k* (*t*) , ∂*U* (*t*)

( ) ( )

= ¶ ¶ ¶+ = ¶ ¶¶ å *m k k k j*

*Ut Ut u t*

Somewhat simplified structure of DHP is shown in Figure 8a. Partial derivatives

( ) ( )

( ) 1

∂*uk* (*t*) are determined from model network, and partial derivatives

Action neural network is trained to minimize Bellman sum, min *J*(*t*) , from the criterion

[ ] ˆ( 1) ( ) [ ] 0 () ()

= × å () () *T A AA t*

[ ] ˆ( 1) ( ) [ ] ( ) () ()

Cross training of critic and action networks results in slower approximation. In order to reduce training time, it is recommended to begin the procedure with action network pre-trained as

¶ +¶ = + ¶ ¶ *A J yt U yt E t*

Globalized Dual heuristic optimal control algorithm uses both algorithms described.

Critic neural network provides approximation of Bellman sum *J*(*t*) and derivative ∂ *J*(*t* + 1)

g

Simplified training procedure of action NN is shown in Figure 8a.

+ =

¶ +¶

¶ ¶ *J yt U yt*

( )

1

*Yt u t y t* (23)

Dual Heuristic Neural Programming Controller for Synchronous Generator

^(*t* + 1) and *λ*(*t* + 1) is determined from model

*yt yt* (24)

*E EtEt* (25)

*yt yt* (26)

∂ *y* ^ *i* (*t* + 1)

http://dx.doi.org/10.5772/58377

∂*u* ^ *k* (*t* + 1) ∂ *y* ^ *j* (*t*)

∂ *y* ^ *j* (*t*) , 117

∂ *y*(*t*)

Model network training is the same as in HDP procedure. The improvement is in the direct calculation of the Bellman sum partial derivative over state space variables ( 1)

(b)

Figure 8. DHP structure, (a) Critic NN training, (b) action NN training **Figure 8.** DHP structure, (a) Critic NN training, (b) action NN training

$$\frac{\partial J\left(t+1\right)}{\partial \hat{\boldsymbol{\gamma}}\_{\boldsymbol{\cdot}}\left(t\right)} = \sum\_{l=1}^{n} \lambda\_{l}\left(t+1\right) \frac{\partial \hat{\boldsymbol{\gamma}}\_{\boldsymbol{\cdot}}\left(t+1\right)}{\partial \hat{\boldsymbol{\gamma}}\_{\boldsymbol{\cdot}}\left(t\right)} + \sum\_{k=1}^{m} \sum\_{l=1}^{n} \lambda\_{l}\left(t+1\right) \frac{\partial \hat{\boldsymbol{\gamma}}\_{\boldsymbol{\cdot}}\left(t+1\right)}{\partial \hat{\boldsymbol{u}}\_{k}\left(t\right)} \frac{\partial \hat{\boldsymbol{u}}\_{k}\left(t+1\right)}{\partial \hat{\boldsymbol{\gamma}}\_{\boldsymbol{\cdot}}\left(t\right)}\tag{22}$$

where

$$
\widehat{\lambda}(t+1) = \frac{\partial \widehat{f\mathbb{L}}\_{\infty}^{\wedge}(t+1)}{\partial \, y(t+1)},
$$

*n* is size of model network output vector *y* ^(*t*), *m* is size of control vector *u* ^(*t*).

Third part of equation (22) is determined from deriving

$$\frac{\partial U\begin{pmatrix} t \\ \end{pmatrix}}{\partial Y\begin{pmatrix} t \end{pmatrix}} = \sum\_{k=1}^{m} \frac{\partial U\begin{pmatrix} t \end{pmatrix}}{\partial u\_k\begin{pmatrix} t \end{pmatrix}} \frac{\partial u\_k\begin{pmatrix} t+1 \\ \end{pmatrix}}{\partial \boldsymbol{\psi}\_f\begin{pmatrix} t \end{pmatrix}} \tag{23}$$

Somewhat simplified structure of DHP is shown in Figure 8a. Partial derivatives ∂ *y* ^ *i* (*t* + 1) ∂ *y* ^ *j* (*t*) , ∂ *y* ^ *i* (*t* + 1) ∂*u* ^ *k* (*t*) , ∂*U* (*t*) ∂*uk* (*t*) are determined from model network, and partial derivatives ∂*u* ^ *k* (*t* + 1) ∂ *y* ^ *j* (*t*) from action network. Prediction of states *y* ^(*t* + 1) and *λ*(*t* + 1) is determined from model network.

Action neural network is trained to minimize Bellman sum, min *J*(*t*) , from the criterion ∂ *J y* ^(*t*) ∂ *y* ^(*t*) ≈0.

Thus, equation (22) can be transformed to

$$\gamma \frac{\partial J\left[\hat{\boldsymbol{\varphi}}(t+1)\right]}{\partial \boldsymbol{\chi}(t)} + \frac{\partial U\left[\boldsymbol{\chi}(t)\right]}{\partial \boldsymbol{\chi}(t)} = 0 \tag{24}$$

Training of Action Network minimizes scalar product:

$$\left\| E\_A \right\| = \sum\_t E\_A^T(t) \cdot E\_A(t) \tag{25}$$

with

( )

4.2 Dual Heuristic Programming

output of critic network.

116 MATLAB Applications for the Practical Engineer

TD

TD

TD

TD

TD

TD

y t( 1) −

u t( )

TD

TD

TD

y t( 1) −

u t( )

TD

y t( )

procedure efficiency by increasing training speed.

Structure of DHP control system is shown in Figure 8.

NETWORK

ACTION

NETWORK

( ) ( ) U t u t ∂ ∂

MODEL y t ˆ( )

u t ˆ( )

MODEL y t ˆ( )

where

∂ *J y* ^(*t* + 1)

∂ *y*

*m* is size of control vector *u*

^(*t* + 1) ,

*n* is size of model network output vector *y*

*λ* ^(*t* + 1)= ( ) ( ) ( )

= = =

l

**Figure 8.** DHP structure, (a) Critic NN training, (b) action NN training

^(*t*).

( ) ( ) ( )

TD

y t( )

( ) ( ) U t u t ∂ ∂

TD

*y t y t ut yt* (22)

( ) 1 1 1 1 ˆ 1 ˆ ˆ 1 1 1 1 ˆ ˆ ˆ ˆ

Figure 8. DHP structure, (a) Critic NN training, (b) action NN training

(b)

Dual Heuristic Programming optimal control algorithm is developed with a goal of increasing

TD

y t ˆ( 1) +

TD

TD

y t ˆ( 1) +

( 1) ( 1) ( 1) () () yt ut t yt yt λ ∂+ ∂+ +⋅ ⋅ ∂ ∂

(a)

( 1) ( 1) ( ) y t t u t λ γ ∂ + +⋅ ⋅ ∂

y t ˆ( )

() () () () Ut ut ut yt ∂ ∂⋅ ∂ ∂

( 1) ( 1) ( ) y t t y t λ ∂ + + ⋅ ∂

TD

TD

TD

Model network training is the same as in HDP procedure. The improvement is in the direct calculation of the Bellman sum partial derivative over state space variables ( 1)

> ¶ + ¶ + ¶ +¶ + =+ + + ¶ ¶ ¶ ¶ å åå *n m n i i k i i j i j k i k j J t y t yt u t t t*

> > ^(*t*),

 l

( )

( )

ACTION

+ +

( ) J t y t ∂ + ∂

CRITIC

+ - - -

CRITIC

CRITIC

as

λ( 1) t +

λ( )t

λ( 1) t +

$$E\_A(t) = \gamma \frac{\partial J\left[\dot{\mathcal{Y}}(t+1)\right]}{\partial \mathbf{y}(t)} + \frac{\partial U\left[\mathbf{y}(t)\right]}{\partial \mathbf{y}(t)}\tag{26}$$

Simplified training procedure of action NN is shown in Figure 8a.

Cross training of critic and action networks results in slower approximation. In order to reduce training time, it is recommended to begin the procedure with action network pre-trained as linear controller.

#### **4.3. Global dual heuristic programming**

Globalized Dual heuristic optimal control algorithm uses both algorithms described.

Critic neural network provides approximation of Bellman sum *J*(*t*) and derivative ∂ *J*(*t* + 1) ∂ *y*(*t*)

Action network is trained using the same procedure as described for HDP and DHP.

point after each disturbance. For those conditions, power system stabilizer (*PSS*) is added to AVR. Based on change of output power and load angle, it modifies voltage regulator setpoint and increases system damping. PSS is configured using linear control theory, and is presently

Dual Heuristic Neural Programming Controller for Synchronous Generator

http://dx.doi.org/10.5772/58377

119

Large disturbances occur after short circuits on the power grid, transmission and generator

Electromechanical oscillations occurring as a consequence of disturbances cause damped oscillations of all electrical and mechanical values (generator voltages and currents, active and reactive power, speed and load angle). Oscillations are also propagated by AVR to excitation

Large system disturbances can cause significant changes to local parameters relevant for generator stability. Due to the linear nature of the voltage regulator system, a change in grid voltage or configuration leads to generator working in non optimal state. Even if the system returns to the same operating conditions, passing through nonlinear states during disturbance causes non optimal control, diminishing the stabilizing effects of the voltage regulator on the system. This can in turn cause another disturbance or, in case of larger units, disruption of the whole power grid. Limitations of the use of linear control systems call for an improved system

There is a lot of research effort directed toward design of excitation control system based on neural networks. However, very few of the demonstrated solutions have been implemented

The goal of the chapter is to develop a practical-to-use system. This is mainly achieved by separating the design on preparation and implementation phase. Almost all algorithm steps are executed in off-line mode, up to the final step of digital voltage regulator implementation.

To implement the voltage regulator described in the chapter following equipment is needed:

**•** automatic voltage regulator with neural network support (no training support is needed)

Described controller is not adaptive. Adaptive critic design (ACD) algorithm is used for

Classic voltage regulators (AVR) equipped with PSS can be very effective in voltage control and power system stabilization if it is used in well configured transmission grid. In such conditions the generator stays in the linear operating range and the controller maintains optimal configuration on change of active and reactive power. Voltage regulator accuracy, balance of reactive power, stabilizing effect during large disturbances and active stabilization

which can provide optimal control of generator excitation system [13]- [15].

a standard part of automatic voltage regulator systems [10].

outages, and connection and disconnection of large loads.

system values.

on real plants.

optimal control.

**•** data acquisition hardware and software

in normal operating conditions are achieved.

**•** personal computer with Matlab and Simulink environment

#### **5. Neural networks based control of excitation voltage**

Main elements of electric power system are generators, transmission and distribution network and consumers. Transmission and distribution network form a power grid, which is usually very large in geographic terms. Grid control is implemented on both local and global level, mostly on the power generation side of the grid. Specific to power grid control is large number of closed control loops governing the power and voltage of synchronous generators. Trans‐ mission network is governed by setting voltage on transformers in order to minimize trans‐ mission losses. Governing of generator power and voltage achieves balance of power as well as reactive power flow on the transmission system. Voltage control also has a large impact on the dynamic system behavior.

Generator power output is governed on the generator unit. For water and steam turbine generators, it defines local power transmission to the grid, while also influencing system frequency. Considering the inertness of the generator systems, power output is not difficult to govern. Voltage control manages the balance between generator and grid voltage, and determines reactive power flow. A change in generator voltage setpoint has impact on system state variables and parameters. Disregarding system transients, a change of generator voltage causes a change on generator busbars, generator current and reactive power.

If we also consider system dynamics, generator voltage change also changes system parame‐ ters influencing system damping. Synchronous generator connected to the power grid is dominantly nonlinear in characteristics. Generator output power is proportional to generator and grid voltage, while the relation to angle between generator and grid voltage is a sinus

function. Maximum theoretical output power is realized with a load angle being *δ* = *π* 2 .

Exceeding of maximum value of load angle causes large disturbances and disconnection from power grid.

Operating point of generator connected to the grid is determined by steady state excitation current. Excitation voltage and current values dictate generator voltage, reactive power and load angle.

Frequent small disturbances expected and part of normal operating conditions, resulting in oscillations of generator electromechanical variables. Especially pronounced is the swinging of generator power, reaching several percents of nominal value, with a frequency around 1Hz. It is caused by low system damping and outer disturbances on the grid. Also present on the grid are oscillations caused by the system characteristic frequencies, usually around 0.1Hz.

Power grid can be represented as a single machine connected to an infinite bus system. For each operating point a characteristic frequency and damping can be determined.

In case of small disturbances with a synchronous generator working around predefined operating point, system can be considered linear, but only if it returns to the same operating point after each disturbance. For those conditions, power system stabilizer (*PSS*) is added to AVR. Based on change of output power and load angle, it modifies voltage regulator setpoint and increases system damping. PSS is configured using linear control theory, and is presently a standard part of automatic voltage regulator systems [10].

Large disturbances occur after short circuits on the power grid, transmission and generator outages, and connection and disconnection of large loads.

Electromechanical oscillations occurring as a consequence of disturbances cause damped oscillations of all electrical and mechanical values (generator voltages and currents, active and reactive power, speed and load angle). Oscillations are also propagated by AVR to excitation system values.

Large system disturbances can cause significant changes to local parameters relevant for generator stability. Due to the linear nature of the voltage regulator system, a change in grid voltage or configuration leads to generator working in non optimal state. Even if the system returns to the same operating conditions, passing through nonlinear states during disturbance causes non optimal control, diminishing the stabilizing effects of the voltage regulator on the system. This can in turn cause another disturbance or, in case of larger units, disruption of the whole power grid. Limitations of the use of linear control systems call for an improved system which can provide optimal control of generator excitation system [13]- [15].

There is a lot of research effort directed toward design of excitation control system based on neural networks. However, very few of the demonstrated solutions have been implemented on real plants.

The goal of the chapter is to develop a practical-to-use system. This is mainly achieved by separating the design on preparation and implementation phase. Almost all algorithm steps are executed in off-line mode, up to the final step of digital voltage regulator implementation.

To implement the voltage regulator described in the chapter following equipment is needed:

**•** data acquisition hardware and software

Action network is trained using the same procedure as described for HDP and DHP.

Main elements of electric power system are generators, transmission and distribution network and consumers. Transmission and distribution network form a power grid, which is usually very large in geographic terms. Grid control is implemented on both local and global level, mostly on the power generation side of the grid. Specific to power grid control is large number of closed control loops governing the power and voltage of synchronous generators. Trans‐ mission network is governed by setting voltage on transformers in order to minimize trans‐ mission losses. Governing of generator power and voltage achieves balance of power as well as reactive power flow on the transmission system. Voltage control also has a large impact on

Generator power output is governed on the generator unit. For water and steam turbine generators, it defines local power transmission to the grid, while also influencing system frequency. Considering the inertness of the generator systems, power output is not difficult to govern. Voltage control manages the balance between generator and grid voltage, and determines reactive power flow. A change in generator voltage setpoint has impact on system state variables and parameters. Disregarding system transients, a change of generator voltage

If we also consider system dynamics, generator voltage change also changes system parame‐ ters influencing system damping. Synchronous generator connected to the power grid is dominantly nonlinear in characteristics. Generator output power is proportional to generator and grid voltage, while the relation to angle between generator and grid voltage is a sinus function. Maximum theoretical output power is realized with a load angle being *δ* = *π*

Exceeding of maximum value of load angle causes large disturbances and disconnection from

Operating point of generator connected to the grid is determined by steady state excitation current. Excitation voltage and current values dictate generator voltage, reactive power and

Frequent small disturbances expected and part of normal operating conditions, resulting in oscillations of generator electromechanical variables. Especially pronounced is the swinging of generator power, reaching several percents of nominal value, with a frequency around 1Hz. It is caused by low system damping and outer disturbances on the grid. Also present on the grid are oscillations caused by the system characteristic frequencies, usually around 0.1Hz. Power grid can be represented as a single machine connected to an infinite bus system. For

In case of small disturbances with a synchronous generator working around predefined operating point, system can be considered linear, but only if it returns to the same operating

each operating point a characteristic frequency and damping can be determined.

2 .

causes a change on generator busbars, generator current and reactive power.

**5. Neural networks based control of excitation voltage**

the dynamic system behavior.

118 MATLAB Applications for the Practical Engineer

power grid.

load angle.

Described controller is not adaptive. Adaptive critic design (ACD) algorithm is used for optimal control.

Classic voltage regulators (AVR) equipped with PSS can be very effective in voltage control and power system stabilization if it is used in well configured transmission grid. In such conditions the generator stays in the linear operating range and the controller maintains optimal configuration on change of active and reactive power. Voltage regulator accuracy, balance of reactive power, stabilizing effect during large disturbances and active stabilization in normal operating conditions are achieved.

If the generator is connected to relatively weak transmission grid, with possible configuration changes during operation, use of neural network based control of excitation voltage leads to better results.

#### **5.1. HDP algorithm implementation**

In order to achieve optimal control of synchronous generator connected to the power grid, it is necessary to select a minimum of two system variables – generator voltage and either output power or angular velocity. Selected variables must reflect system oscillations during distur‐ bances. It is also necessary to choose cost function that is minimized during training. System optimization is achieved by minimization of squared generator voltage governing error *ΔU* (*t*)=*UG REF* −*UG* and minimization of output power change *ΔP* obtained from:

$$\begin{split} U(t) &= \left( K\_{\text{ul}} \cdot \Delta U\_G(t) + K\_{\text{ul}2} \cdot \Delta U\_G(t-\text{l}) + K\_{\text{ul}3} \cdot \Delta U\_G(t-\text{2}) \right)^2 + \\ & \left( K\_{P1} \cdot \Delta P(t) + K\_{P1} \cdot \Delta P(t-\text{l}) + K\_{P1} \cdot \Delta P(t-\text{2}) \right)^2 \end{split} \tag{27}$$

1/z 1/z

. It was reduced during training to *η* =102

.

RNN Model

[Ug] [Pg]

Dual Heuristic Neural Programming Controller for Synchronous Generator

[Ug] [Pg]

A third order nonlinear mathematical model was used as a synchronous generator [20]. Synchronous generator is modeled as a single machine connected to an infinite bus. Although the model ignores initial magnetic flows, it proves to be a good approximation of an actual

Selected mathematical model, used with well conditioned state space variables and described methods, achieves good numerical stability during simulation, allowing for relatively large sample time (TD=0.02s). Described model of synchronous generator is implemented in C-

For model identification, recurrent neural network is used, with five neurons in hidden layer and two neurons in output layer. Hidden layer is implemented with full recursion and *tansig*

Training was completed in ~10000 steps. Results are shown in Figure 10. Signals *UG* NN(0) and *UG* NN(1) represent generator voltage at the start and at the end of RNN training. Figure 10a compares responses of control error of generator model and neural network. Voltage setpoint is generated as a pulse signal with a magnitude of ± 10%. In figure 10b, response to the change of generator power is shown. Figure 10c shows squared sum of RNN error during

In the same experiment efficiency of the algorithm was verified. Training was stopped after 100000 steps. For the rest of the simulation, network output is calculated without the change

generator. A predictor-corrector method is used to solve differential equations.

[Ug] [Pg] Evaulation

http://dx.doi.org/10.5772/58377

121

Model of Synchronous generator tm TD

**Figure 9.** Model identification block diagram

language S-function and used for further simulations.

Initial learning rate value *η* (21) is set to *η* =105

of weight coefficient, thus eliminating the recency effect.

activation function. Output layer uses linear activation function.

[Pg]

[Ug] [Pg]

[Ug]

V

training.

u(t)

It is important to notice that the two demands are in conflict. Increase in gain in main control loop reduces static control error, but also reduces system damping and increases power output deviation *ΔP*. Use of cost function (27) ensures training process that increases system damping without sacrificing controller responsiveness.

System identification is carried by training model neural network. Figure 8 shows training process of model neural network. Identification of nonlinear system is made by series-parallel structure.

For ACD algorithm optimal control two output state space variables are chosen – generator voltage and output power, *y* = *ΔUG*;*ΔP* . Vector *y* is used to form the cost function. Input vector is

$$\mathbf{u}(t) = \begin{bmatrix} \mathbf{y}\_r(t) \ \mathbf{y}\_r(t-l) \ \mathbf{y}\_r(t-2) \ \mathbf{y}(t-l) \ \mathbf{y}(t-2) \ \mathbf{y}(t-3) \end{bmatrix}^T \tag{28}$$

where

*yr* =*UGREF* represents generator voltage setpoint,

*uf* – exciter control signal,

*y* = *ΔUG ΔP T* is system output.

Model NN has 9 inputs and 2 outputs.

Figure 9 shows simulation block diagram of identification procedure, implemented in Matlab and Simulink environment.

**Figure 9.** Model identification block diagram

If the generator is connected to relatively weak transmission grid, with possible configuration changes during operation, use of neural network based control of excitation voltage leads to

In order to achieve optimal control of synchronous generator connected to the power grid, it is necessary to select a minimum of two system variables – generator voltage and either output power or angular velocity. Selected variables must reflect system oscillations during distur‐ bances. It is also necessary to choose cost function that is minimized during training. System optimization is achieved by minimization of squared generator voltage governing error

*ΔU* (*t*)=*UG REF* −*UG* and minimization of output power change *ΔP* obtained from:

12 3

*Ut K U t K U t K U t*

( ) ( ) ( 1) ( 2) ( ) ( 1) ( 2)

*uG u G uG*

( )

×D + ×D - + ×D -

11 1

*K Pt K Pt K Pt*

*PP P*

without sacrificing controller responsiveness.

*yr* =*UGREF* represents generator voltage setpoint,

( )

= ×D + ×D - + ×D - +

It is important to notice that the two demands are in conflict. Increase in gain in main control loop reduces static control error, but also reduces system damping and increases power output deviation *ΔP*. Use of cost function (27) ensures training process that increases system damping

System identification is carried by training model neural network. Figure 8 shows training process of model neural network. Identification of nonlinear system is made by series-parallel

For ACD algorithm optimal control two output state space variables are chosen – generator voltage and output power, *y* = *ΔUG*;*ΔP* . Vector *y* is used to form the cost function. Input

Figure 9 shows simulation block diagram of identification procedure, implemented in Matlab

( ) ( ) ( 1) ( 2) ( 1) ( 2) ( 3) = - - -- - [ ]

2

*T*

(27)

2

*rr r ut y t y t y t yt yt yt* (28)

better results.

structure.

vector is

where

*uf* – exciter control signal,

and Simulink environment.

*y* = *ΔUG ΔP T* is system output.

Model NN has 9 inputs and 2 outputs.

**5.1. HDP algorithm implementation**

120 MATLAB Applications for the Practical Engineer

A third order nonlinear mathematical model was used as a synchronous generator [20]. Synchronous generator is modeled as a single machine connected to an infinite bus. Although the model ignores initial magnetic flows, it proves to be a good approximation of an actual generator. A predictor-corrector method is used to solve differential equations.

Selected mathematical model, used with well conditioned state space variables and described methods, achieves good numerical stability during simulation, allowing for relatively large sample time (TD=0.02s). Described model of synchronous generator is implemented in Clanguage S-function and used for further simulations.

For model identification, recurrent neural network is used, with five neurons in hidden layer and two neurons in output layer. Hidden layer is implemented with full recursion and *tansig* activation function. Output layer uses linear activation function.

Initial learning rate value *η* (21) is set to *η* =105 . It was reduced during training to *η* =102 . Training was completed in ~10000 steps. Results are shown in Figure 10. Signals *UG* NN(0) and *UG* NN(1) represent generator voltage at the start and at the end of RNN training. Figure 10a compares responses of control error of generator model and neural network. Voltage setpoint is generated as a pulse signal with a magnitude of ± 10%. In figure 10b, response to the change of generator power is shown. Figure 10c shows squared sum of RNN error during training.

In the same experiment efficiency of the algorithm was verified. Training was stopped after 100000 steps. For the rest of the simulation, network output is calculated without the change of weight coefficient, thus eliminating the recency effect.

Simulation uses expressions (21), (22) and (23).

to system stability during large disturbances.

output signals and partial derivatives of outputs ( ∂ *yi*

language.

trained in parallel.

> 0.4 0.5 0.6 0.7

d[rad]

power; c) load angle

DP [p.u.]

DUG[p.u.]

Neural network training is implemented in Matlab S-functions, written in C programming

Input of the S function is formed of system signals and desired value. S-function returns RNN

Training results are shown in Figure 12. Figure 12a compares generator voltage response using conventional automatic voltage regulator (AVR) to neural network based controller trained using HDP algorithm. It can be seen that voltage rate of change was not degraded as a result of system damping. It is important to emphasize that generator voltage rate of change is key

Figure 12b compares active power output oscillation of conventional and HDP controller.

(a)

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

t, s

**Figure 12.** Response comparison of conventional AVR and HDP controller. a) generator voltage; b) generator active

t, s

(c)

t, s

(b)

Figure 12c shows load angle, which is indicative of system stability (*δ* <

is significantly increased by use of neural network controller.

(*t*) ∂*uj*

Dual Heuristic Neural Programming Controller for Synchronous Generator

(*t*) ). Critic and action networks are

http://dx.doi.org/10.5772/58377

123

*π*

2 ). System damping

RNN AVR D UGREF

> RNN AVR

AVR RNN

Figure 10.Training results of model neural network **Figure 10.** Training results of model neural network

Once the model neural network is trained, it is possible to create simulation for training critic and action neural networks. Simulink model of the simulation is shown in Figure 11. Once the model neural network is trained, it is possible to create simulation for training critic and action neural networks. Simulink model of the simulation is shown in Figure 11.

Figure 11.Simulink model for training action and critic networks **Figure 11.** Simulink model for training action and critic networks

Simulation uses expressions (21), (22) and (23).

Figure 10.Training results of model neural network

k

1 2 3 4 5 6 7 8 9 10 11

Learn

[Ug]

[Pg]

[dY2du]

[dY1du]

T D

[dY2du] [dY1du]

Action learn

User criteria U(t)

[dUdu]

[dY2du] [dY1du]

NN-Critic HDP

x( t) x(t+ 1) x 104

[dUdu]

P NN(0) P P NN(1)

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

t, s

t, s

(b)

(c)

(a)

∆ UGRE:

122 MATLAB Applications for the Practical Engineer

F ∆UG NN(0) ∆ UG ∆ UG NN(1)

0 0.002 0.004 0.006 0.008 0.01

**Figure 10.** Training results of model neural network

[Ug]

[Pg]

SG

Pg Ug th RNiz

e tm Ru T D

SG Referent

> Pg Ug th

e tm

[Ug]

[Pg]

∆P [p.u.]

Eror

UG REF

Pref

∆UG [p.u.]

Once the model neural network is trained, it is possible to create simulation for training critic

T D

Once the model neural network is trained, it is possible to create simulation for training critic

Figure 11.Simulink model for training action and critic networks

Learn

T D

NN Model

NN action

**Figure 11.** Simulink model for training action and critic networks

Ug P

and action neural networks. Simulink model of the simulation is shown in Figure 11.

and action neural networks. Simulink model of the simulation is shown in Figure 11.

Neural network training is implemented in Matlab S-functions, written in C programming language.

Input of the S function is formed of system signals and desired value. S-function returns RNN output signals and partial derivatives of outputs ( ∂ *yi* (*t*) ∂*uj* (*t*) ). Critic and action networks are trained in parallel.

Training results are shown in Figure 12. Figure 12a compares generator voltage response using conventional automatic voltage regulator (AVR) to neural network based controller trained using HDP algorithm. It can be seen that voltage rate of change was not degraded as a result of system damping. It is important to emphasize that generator voltage rate of change is key to system stability during large disturbances.

Figure 12b compares active power output oscillation of conventional and HDP controller. Figure 12c shows load angle, which is indicative of system stability (*δ* < *π* 2 ). System damping is significantly increased by use of neural network controller.

**Figure 12.** Response comparison of conventional AVR and HDP controller. a) generator voltage; b) generator active power; c) load angle

HDP algorithm is also verified on voltage control of synchronous generator from Simulink SimPower System block library (Synchronous Machine).

Figure 13b shows comparison of voltage transient in the condition generator voltage reduced to 90% of nominal value. Rate of change of generator voltage is maintained even during such

Dual Heuristic Neural Programming Controller for Synchronous Generator

http://dx.doi.org/10.5772/58377

125

Figure 13d shows conventional automatic voltage regulator (AVR) operating on the edge of stability, while the DHP controller is completely stable. Comparison of load angle values demonstrates significantly better behavior of DHP controller. On large load angles DHP controller maintains stability while conventional AVR provides low system damping, leading

Neural network HDP voltage controllers were developed using DHP algorithm. Experimental results showed that DHP algorithm is more efficient than ADC algorithm. Implementation of

Application of developed S-functions enables simple use of described algorithms, as well as

DHP optimal control algorithm has also been implemented using Matlab Real Time Windows Target Toolbox (RTWT), based on [18] and [19]. Simulation hardware consists of desktop PC

Use of RTWT platform enables direct implementation of developed controllers on the real

In Figure 14, a system consisting of synchronous generator and RTWT controller is shown. Signal acquisition is achieved using DAQ 6062E Input, with AD conversion of the input values. Line voltage and phase current signals are processed in "Clark" S-function, using Clarke transformation. Transformed values of generator voltage *UG* and generator power *P* are used as control feedbacks. Digital to analog conversion is performed by output module DAQ 6062 Output. All simulation blocks are part of a standard Simulink library except Clark transfor‐

ACD algorithm is implemented in two steps. Classic voltage controller is used in the first step

During testing, generator active power and voltage were changed for a wide range of values. A pulse signal was superimposed to the generator voltage set-point value, resulting in large disturbances on the generator. To obtain measurement records RTWT Signal Visualization was used, with measured data saved in Workspace Simulink block *yP*. Recorded data was

A RNN was used with one hidden layer containing 5 tansig neurons. Model NN had 9 inputs

equipped with National Instruments data acquisition card NI DAQ 6024.

(Figure 14). In the second step it is replaced with NN controller.

and 2 outputs. Results of NN training are displayed in Figure 15.

extreme operating conditions.

to instability and possible outages.

system.

mation and RNN S-functions.

used to train neural network in offline mode.

GDHP does not bring significant improvements over DHP.

experimenting with different function parameters.

**6. Real time windows target implementation**

Model network is designed with one hidden recursion layer with 5 neurons. Critic and action networks use the same structure, but with four neurons in hidden layer.

HDP algorithm shown in Figure 11 is applied. Training of neural networks is performed in two steps– in first step action NN is trained to mimic classic voltage regulator, and in second step DHP algorithm is used to train both action and critic NN. It is important to get initial values of action NN in first step, as DHP needs network partial derivatives in training critic NN. Figure 13 shows training results.

Figure 13a shows comparison of classic and NN controller responses. Generator voltage rate of change is almost the same for both controllers.

**Figure 13.** Response comparison of conventional AVR and HDP controller a), b) generator voltage; c), d) generator active power; e), f) load angle

Figure 13c and 13e show comparison of output power and load angle change. It is obvious that the HDP neural network controller shows significant improvement over classic regulator.

Figure 13b shows comparison of voltage transient in the condition generator voltage reduced to 90% of nominal value. Rate of change of generator voltage is maintained even during such extreme operating conditions.

Figure 13d shows conventional automatic voltage regulator (AVR) operating on the edge of stability, while the DHP controller is completely stable. Comparison of load angle values demonstrates significantly better behavior of DHP controller. On large load angles DHP controller maintains stability while conventional AVR provides low system damping, leading to instability and possible outages.

Neural network HDP voltage controllers were developed using DHP algorithm. Experimental results showed that DHP algorithm is more efficient than ADC algorithm. Implementation of GDHP does not bring significant improvements over DHP.

Application of developed S-functions enables simple use of described algorithms, as well as experimenting with different function parameters.

### **6. Real time windows target implementation**

HDP algorithm is also verified on voltage control of synchronous generator from Simulink

Model network is designed with one hidden recursion layer with 5 neurons. Critic and action

HDP algorithm shown in Figure 11 is applied. Training of neural networks is performed in two steps– in first step action NN is trained to mimic classic voltage regulator, and in second step DHP algorithm is used to train both action and critic NN. It is important to get initial values of action NN in first step, as DHP needs network partial derivatives in training critic

Figure 13a shows comparison of classic and NN controller responses. Generator voltage rate

> -0.2 -0.1 0 0.1 0.2

0 0.2 0.4 0.6 0.8 1 1.2

∆ δ [rad]

**Figure 13.** Response comparison of conventional AVR and HDP controller a), b) generator voltage; c), d) generator

Figure 13c and 13e show comparison of output power and load angle change. It is obvious that the HDP neural network controller shows significant improvement over classic regulator.

(e) (f)

∆UGREF NN AVR

∆ P [p.u.]

(c) (d)

∆ UG [p.u.]

∆ UGREF ∆ UG NN ∆ UG AVR

0 1 2 3 4 5 6 7 8 9 10

(b)

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

t, s

t, s

∆UGREF NN P AVR

t, s

∆ UGREF ∆ UG NN ∆ UG AVR

> ∆UGREF NN AVR

SimPower System block library (Synchronous Machine).

NN. Figure 13 shows training results.

124 MATLAB Applications for the Practical Engineer

> 0 0.2 0.4 0.6 0.8

active power; e), f) load angle

∆ δ [rad]

∆ P [p.u.]

∆ VG [p.u.]

of change is almost the same for both controllers.

(a)

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

t, s

t, s

∆UGREF NN AVR

t, s

networks use the same structure, but with four neurons in hidden layer.

DHP optimal control algorithm has also been implemented using Matlab Real Time Windows Target Toolbox (RTWT), based on [18] and [19]. Simulation hardware consists of desktop PC equipped with National Instruments data acquisition card NI DAQ 6024.

Use of RTWT platform enables direct implementation of developed controllers on the real system.

In Figure 14, a system consisting of synchronous generator and RTWT controller is shown. Signal acquisition is achieved using DAQ 6062E Input, with AD conversion of the input values. Line voltage and phase current signals are processed in "Clark" S-function, using Clarke transformation. Transformed values of generator voltage *UG* and generator power *P* are used as control feedbacks. Digital to analog conversion is performed by output module DAQ 6062 Output. All simulation blocks are part of a standard Simulink library except Clark transfor‐ mation and RNN S-functions.

ACD algorithm is implemented in two steps. Classic voltage controller is used in the first step (Figure 14). In the second step it is replaced with NN controller.

During testing, generator active power and voltage were changed for a wide range of values. A pulse signal was superimposed to the generator voltage set-point value, resulting in large disturbances on the generator. To obtain measurement records RTWT Signal Visualization was used, with measured data saved in Workspace Simulink block *yP*. Recorded data was used to train neural network in offline mode.

A RNN was used with one hidden layer containing 5 tansig neurons. Model NN had 9 inputs and 2 outputs. Results of NN training are displayed in Figure 15.

system.

6 Real Time Windows Target Implementation

Figure 14.Automatic voltage control of synchronous generator using RTWT platform

A pulse signal was superimposed to the generator voltage set-point value, resulting in large disturbances on the generator. To obtain measurement records RTWT Signal Visualization

DHP optimal control algorithm has also been implemented using Matlab Real Time Windows

In Figure 14, a system consisting of synchronous generator and RTWT controller is shown. Signal acquisition is achieved using DAQ 6062E Input, with AD conversion of the input values. Line voltage and phase current signals are processed in "Clark" S-function, using Clarke transformation. Transformed values of generator voltage UG and generator power P are used as control feedbacks. Digital to analog conversion is performed by output module DAQ 6062 Output. All simulation blocks are part of a standard Simulink library except Clark

> Obtained model NN was used to train critic and action RNN using procedure described in 3. The trained action NN was transferred to Simulink block (Figure 14) and prepared for use in real time. Using RTWT module system was run in real-time mode, with trained action RNN

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10 20

t, s

0 1 2 3 4 5 6 7 8 9 10

**Figure 16.** Classic controller and HDP NN controller – system responses: a) generator voltage; b) generator active pow‐

Simulation results are shown in Figure 16, comparing classic and DHP voltage controllers. Figure 16a represents comparison between generator voltage responses and shows satisfactory rate of generator voltage change and high control accuracy. In Figure 16b, generator active power outputs are compared. Considerable damping was achieved using DHP controller. Improved system stability during large disturbances is visible in Figure 16c, where the

RTWT controller was implemented in sample-time mode with sample time set to

Matlab Simulink RTWT Toolbox simplifies the transfer of developed procedures to the

maximum load angle deviation is significantly lower for DHP controller.

*s*. Control loop sample time is *TD* =0.02*s*.

t, s

t, s

(b)

(a)

Dual Heuristic Neural Programming Controller for Synchronous Generator

(c)

UGREF UG AVR UG NN

127

http://dx.doi.org/10.5772/58377

P AVR P NN

> δ AVR δ NN

as a voltage controller.

0.6 0.65 0.7 0.75 0.8 0.85 0.9

er; c) load angle

*Ts* =2.5⋅10−4

hardware platform.

δ [ ° ]

P [kW]

UG [p.u.]

Target Toolbox (RTWT), based on [18] and [19]. Simulation hardware consists of desktop PC equipped with National Instruments data acquisition card NI DAQ 6024. Use of RTWT platform enables direct implementation of developed controllers on the real

step (Figure 14). In the second step it is replaced with NN controller.

During testing, generator active power and voltage were changed for a wide range of values. **Figure 14.** Automatic voltage control of synchronous generator using RTWT platform

**Figure 15.** Results of Model NN training a) generator voltage; b) generator active power

Obtained model NN was used to train critic and action RNN using procedure described in 3. The trained action NN was transferred to Simulink block (Figure 14) and prepared for use in real time. Using RTWT module system was run in real-time mode, with trained action RNN as a voltage controller.

**Figure 16.** Classic controller and HDP NN controller – system responses: a) generator voltage; b) generator active pow‐ er; c) load angle

0 2 4 6 8 10 12 14 16 18 20 -0.3

(b)

(a)

A RNN was used with one hidden layer containing 5 tansig neurons. Model NN had 9 inputs

Figure 14.Automatic voltage control of synchronous generator using RTWT platform

During testing, generator active power and voltage were changed for a wide range of values. A pulse signal was superimposed to the generator voltage set-point value, resulting in large disturbances on the generator. To obtain measurement records RTWT Signal Visualization was used, with measured data saved in Workspace Simulink block yP. Recorded data was

6 Real Time Windows Target Implementation

transformation and RNN S-functions.

126 MATLAB Applications for the Practical Engineer

system.

DHP optimal control algorithm has also been implemented using Matlab Real Time Windows

In Figure 14, a system consisting of synchronous generator and RTWT controller is shown. Signal acquisition is achieved using DAQ 6062E Input, with AD conversion of the input values. Line voltage and phase current signals are processed in "Clark" S-function, using Clarke transformation. Transformed values of generator voltage UG and generator power P are used as control feedbacks. Digital to analog conversion is performed by output module DAQ 6062 Output. All simulation blocks are part of a standard Simulink library except Clark

ACD algorithm is implemented in two steps. Classic voltage controller is used in the first

step (Figure 14). In the second step it is replaced with NN controller.

Target Toolbox (RTWT), based on [18] and [19]. Simulation hardware consists of desktop PC equipped with National Instruments data acquisition card NI DAQ 6024. Use of RTWT platform enables direct implementation of developed controllers on the real

t, s

0 2 4 6 8 10 12 14 16 18 20

DUGREF DP DP NN

DUGREF DUG DUG NN

t, s

**Figure 15.** Results of Model NN training a) generator voltage; b) generator active power

used to train neural network in offline mode.

and 2 outputs. Resultsof NN training are displayed in Figure 15.

**Figure 14.** Automatic voltage control of synchronous generator using RTWT platform

DPDUG[p.u.]

DUG[p.u.]

> Simulation results are shown in Figure 16, comparing classic and DHP voltage controllers. Figure 16a represents comparison between generator voltage responses and shows satisfactory rate of generator voltage change and high control accuracy. In Figure 16b, generator active power outputs are compared. Considerable damping was achieved using DHP controller. Improved system stability during large disturbances is visible in Figure 16c, where the maximum load angle deviation is significantly lower for DHP controller.

> RTWT controller was implemented in sample-time mode with sample time set to *Ts* =2.5⋅10−4 *s*. Control loop sample time is *TD* =0.02*s*.

> Matlab Simulink RTWT Toolbox simplifies the transfer of developed procedures to the hardware platform.

#### **7. Conclusion**

In the chapter, a complete process of development and implementation of neural network based controller is described. Optimal control is achieved by training neural network using data obtained from real system. Developed algorithms were tested on voltage control of synchronous generator connected to the power grid.

end

S-function for RNN training with one hidden layer.

input. Training is performed using Kalman filter.

#define S\_FUNCTION\_LEVEL 2 #include "simstruc.h" #include <math.h> #include <string.h>

#define TS 1 //sample time #define N\_ULD 3 //number of inputs

#define NUL (N\_ULD\*N\_TD) #define NIZ 1 //number outputs

#define W2\_j (W1+1) #define W1\_ij (NUL+W1+1)\*W1 #define W2\_ij NIZ\*W2\_j #define N\_W1 0

#define N\_W2 N\_W1+W1\_ij #define SIZE\_P (W1\_ij+W2\_ij) #define N\_P N\_W2+W2\_ij

#define N\_dadWR NH\_1+W1

#define W1\_j (NUL+W1+1) //

#define NH\_1 N\_P+ SIZE\_P\*SIZE\_P

#define N\_OUT N\_dadWR+W1\*W1\_ij #define N\_dnet N\_OUT+NIZ

end,

%end function dadWR

vision: 1.12 \$ \*/

end,

end,

**Appendix B**

dadWR(l,(i-1)\*(Ni+NW1)+j)=Ia(i)\*((l==i)\*p(j)+dadar\_temp);

The first part defines sample time, number of inputs, outputs and neurons. Input is formed of values used for training. Output 1 is RNN output. Output 2 is partial derivative of output over

puts \* Mato Miskovic Imotica \*/#define S\_FUNCTION\_NAME model

#define W1 5 // number of neurons in the hidden layer

#define N\_TD 1 //number of steps backwards (t-1), (t-3),(t-3), (t-N\_TD)

/\* File : model.c \* Abstract: \* For more details about S-functions, see simulink/src/sfuntmpl\_doc.c. \* Copyright 1990-2002 The MathWorks, Inc. \* \$Re-

function is used for learning recurrent neural network (RNN) in programming environment Matlab Simulink \* Learning RNN is implemented using the Kalman filter \* This program supports RNN with one hidden layer \* In the first part defines the time discretization, number of inputs, number of outputs and \* the number of neurons in a single hidden layer \* Simultaneously with learning RNN this function account the partial derivative of output per inputs \* Port number of the s-function is equal to the number of inputs in RNN + number of out-

/\* ================model===================== \* This s-

Dual Heuristic Neural Programming Controller for Synchronous Generator

http://dx.doi.org/10.5772/58377

129

Optimal control algorithm is implemented in Matlab and Simulink environment. Neural network training is implemented in C programming language Matlab S-function. S-function execution time was sufficiently low to achieve real time performance.

Developed algorithm is tested using Matlab Real Time Windows Target Toolbox, using the same models with minimal modifications.

Experiments performed on laboratory system show possibility of building optimal controller based on special network training procedure with no need for adaptive features.

Obtained simulation results show advantages of neural networks based controller over classic controllers. Developed procedure can relatively easily be adapted for use on real systems.

#### **Appendices**

#### **Appendix A**

Matlab function for calculating hidden layer derivatives of outputs over weight coefficients of hidden layer.

```
function[dadWR,a]=dadWR(WW1,WW2,p,a_1,dadWR_1,NW1,NW2);
%function determines the partial derivative output recursive hidden%layers to
weight coefficients matrix WW1% WW1 matrix of hidden layer% WW2 matrix of out-
put layer% NW1 number of neurons in the hidden layer% NW2 length of the output
vector neural networks% Ni length of the input vector neural networks% J tar-
get output value of the neural network% a output value of hidden layer neurons
% p input vector of the hidden layer of recurrent neural network
Ni=length(p)-NW1;
a_temp=WW1*[p; 1];
a=tansig(a_temp);
Ia=ones(NW1,1)-(a.*a); %value of derivative tansig activation function of neu-
rons -
%y=WW2*a;
%J=J-y;
for l=1:NW1
   for i=1:NW1
       for j=1:Ni+NW1
 dadar_temp=0;
           for k=1:NW1;
 dadar_temp=dadar_temp+WW1(l,Ni+k)*dadWR_1(k,(i-1)*Ni+j);
```

```
end
 dadWR(l,(i-1)*(Ni+NW1)+j)=Ia(i)*((l==i)*p(j)+dadar_temp);
       end,
   end,
end,
%end function dadWR
```
#### **Appendix B**

**7. Conclusion**

128 MATLAB Applications for the Practical Engineer

**Appendices**

**Appendix A**

hidden layer.

synchronous generator connected to the power grid.

same models with minimal modifications.

Ni=length(p)-NW1; a\_temp=WW1\*[p; 1]; a=tansig(a\_temp);

for i=1:NW1

for j=1:Ni+NW1 dadar\_temp=0;

for k=1:NW1;

rons - %y=WW2\*a; %J=J-y; for l=1:NW1

execution time was sufficiently low to achieve real time performance.

In the chapter, a complete process of development and implementation of neural network based controller is described. Optimal control is achieved by training neural network using data obtained from real system. Developed algorithms were tested on voltage control of

Optimal control algorithm is implemented in Matlab and Simulink environment. Neural network training is implemented in C programming language Matlab S-function. S-function

Developed algorithm is tested using Matlab Real Time Windows Target Toolbox, using the

Experiments performed on laboratory system show possibility of building optimal controller

Obtained simulation results show advantages of neural networks based controller over classic controllers. Developed procedure can relatively easily be adapted for use on real systems.

Matlab function for calculating hidden layer derivatives of outputs over weight coefficients of

% p input vector of the hidden layer of recurrent neural network

%function determines the partial derivative output recursive hidden%layers to weight coefficients matrix WW1% WW1 matrix of hidden layer% WW2 matrix of output layer% NW1 number of neurons in the hidden layer% NW2 length of the output vector neural networks% Ni length of the input vector neural networks% J target output value of the neural network% a output value of hidden layer neurons

Ia=ones(NW1,1)-(a.\*a); %value of derivative tansig activation function of neu-

dadar\_temp=dadar\_temp+WW1(l,Ni+k)\*dadWR\_1(k,(i-1)\*Ni+j);

based on special network training procedure with no need for adaptive features.

function[dadWR,a]=dadWR(WW1,WW2,p,a\_1,dadWR\_1,NW1,NW2);

S-function for RNN training with one hidden layer.

The first part defines sample time, number of inputs, outputs and neurons. Input is formed of values used for training. Output 1 is RNN output. Output 2 is partial derivative of output over input. Training is performed using Kalman filter.

```
/* File : model.c * Abstract: * For more details about S-functions, see
simulink/src/sfuntmpl_doc.c. * Copyright 1990-2002 The MathWorks, Inc. * $Re-
vision: 1.12 $ */
                       /* ================model===================== * This s-
function is used for learning recurrent neural network (RNN) in programming en-
vironment Matlab Simulink * Learning RNN is implemented using the Kalman
filter * This program supports RNN with one hidden layer * In the first part
defines the time discretization, number of inputs, number of outputs and * the
number of neurons in a single hidden layer * Simultaneously with learning RNN
this function account the partial derivative of output per inputs * Port num-
ber of the s-function is equal to the number of inputs in RNN + number of out-
puts * Mato Miskovic Imotica */#define S_FUNCTION_NAME model
#define S_FUNCTION_LEVEL 2
#include "simstruc.h"
#include <math.h>
#include <string.h>
#define TS 1 //sample time
#define N_ULD 3 //number of inputs
#define N_TD 1 //number of steps backwards (t-1), (t-3),(t-3), (t-N_TD)
#define NUL (N_ULD*N_TD)
#define NIZ 1 //number outputs
#define W1 5 // number of neurons in the hidden layer
#define W1_j (NUL+W1+1) //
#define W2_j (W1+1)
#define W1_ij (NUL+W1+1)*W1
#define W2_ij NIZ*W2_j
#define N_W1 0
#define N_W2 N_W1+W1_ij
#define SIZE_P (W1_ij+W2_ij)
#define N_P N_W2+W2_ij
#define NH_1 N_P+ SIZE_P*SIZE_P
#define N_dadWR NH_1+W1
#define N_OUT N_dadWR+W1*W1_ij
#define N_dnet N_OUT+NIZ
```

```
#define N_ON N_dnet+NIZ*N_ULD
#define NX N_ON+10 //
#define ETTAMI 10.0 //final value of the learning rate
#define ETTAQ .051 // Q
#define ITER_UCI 20000 //number of steps to stop learning RNN
#define NDISCSTATES NX
static
                  void ttsig(real_T x, real_T *y) //tansig
{
 *y=2.0/(1.0+exp(-2.0*x))-1.0;
}
//AxB=C multiplying matrices in line form
static
                    void aputab_c(real_T *aa, real_T *bb, real_T *cc, int_T
a_r, int_T a_c, int_T b_c){
 int_T i, j, k;
 real_T priv;
   for(i=0; i<a_r ;i++){
       for(j=0;j<b_c;j++){
 priv=0;
          for(k=0;k<a_c;k++) {
 priv=priv+aa[k+i*a_c]*bb[j+b_c*k];
 }
 cc[i*b_c+j]=priv;
 }
 }
}
//XxB=C multiplying matrices in line form, x from SimStruc
static
                    void xputab_c_pok(SimStruct *S, real_T *bb, real_T *cc,
int_T a_r, int_T a_c, int_T b_c, int_T x_p){
 int_T i, j, k;
 real_T priv;
 real_T *x = ssGetRealDiscStates(S);
   for(i=0; i<a_r ;i++){
       for(j=0;j<b_c;j++){
 priv=0;
          for(k=0;k<a_c;k++) {
 priv=priv+x[x_p+k+i*a_c]*bb[j+b_c*k];
 }
 cc[i*b_c+j]=priv;
 }
 }
}
//AxX=C multiplying matrices in line form, x from SimStruc
static
                    void aputax_c_pok(SimStruct *S, real_T *aa, real_T *cc,
int_T a_r, int_T a_c, int_T b_c, int_T x_p){
 int_T i, j, k;
 real_T priv;
 real_T *x = ssGetRealDiscStates(S);
```
for(i=0; i<a\_r ;i++){ for(j=0;j<b\_c;j++){

cc[i\*b\_c+j]=priv;

for(i=0; i<a\_r ;i++){ for(j=0;j<b\_c;j++){

x[x\_p+i\*b\_c+j]=priv;

for(i=0; i<a\_c ;i++){

for(l=0;l<n;l++){ p\_old2=l\*n+l;

p\_old=m\*n+l;

 aa[p\_old2]=1/(aa[p\_old2]); for(m=0;m<n;m++) if(m!=l){

for(m=0;m<n;m++){

aa[p\_old]=aa[p\_old]\*aa[p\_old2];

for(p=0;p<n;p++){ if(m!=l){

if(p!=l){

for(j=0; j<a\_r ;j++){ bb[i\*a\_r+j]=aa[i+j\*(a\_c)]; }}

for(k=0;k<a\_c;k++) {

int\_T a\_r, int\_T a\_c, int\_T b\_c, int\_T x\_p){

real\_T \*x = ssGetRealDiscStates(S);

for(k=0;k<a\_c;k++) { priv=priv+aa[k+i\*a\_c]\*bb[j+b\_c\*k];

priv=priv+aa[k+i\*a\_c]\*x[x\_p+j+b\_c\*k]; }

//AxB=X multiplying matrices in line form, x to SimStruc

void aputab\_x\_pok(SimStruct \*S, real\_T \*aa, real\_T \*bb,

Dual Heuristic Neural Programming Controller for Synchronous Generator

http://dx.doi.org/10.5772/58377

131

void transp(real\_T \*aa, real\_T \*bb, int\_T a\_r, int\_T a\_c)

void invmatrix(real\_T \*aa, int\_T n){ // inverse matrixin

priv=0;

 int\_T i, j, k; real\_T priv;

priv=0;

{ //transposed matrix int\_T i, j;

}

 } } } static

} static

inline form int\_T l, m, p; int\_T pos, p\_old; int\_T p\_old2;

}

 } } }

static

```
for(i=0; i<a_r ;i++){
      for(j=0;j<b_c;j++){
 priv=0;
          for(k=0;k<a_c;k++) {
 priv=priv+aa[k+i*a_c]*x[x_p+j+b_c*k]; }
 cc[i*b_c+j]=priv;
 }
 }
}
//AxB=X multiplying matrices in line form, x to SimStruc
static
                  void aputab_x_pok(SimStruct *S, real_T *aa, real_T *bb,
int_T a_r, int_T a_c, int_T b_c, int_T x_p){
 int_T i, j, k;
 real_T priv;
 real_T *x = ssGetRealDiscStates(S);
   for(i=0; i<a_r ;i++){
      for(j=0;j<b_c;j++){
 priv=0;
          for(k=0;k<a_c;k++) {
 priv=priv+aa[k+i*a_c]*bb[j+b_c*k];
 }
 x[x_p+i*b_c+j]=priv;
 }
 }
}
static
                  void transp(real_T *aa, real_T *bb, int_T a_r, int_T a_c)
{ //transposed matrix
 int_T i, j;
   for(i=0; i<a_c ;i++){
      for(j=0; j<a_r ;j++){
 bb[i*a_r+j]=aa[i+j*(a_c)]; }}
}
static
                  void invmatrix(real_T *aa, int_T n){ // inverse matrixin
inline form
 int_T l, m, p;
 int_T pos, p_old;
 int_T p_old2;
   for(l=0;l<n;l++){
 p_old2=l*n+l;
 aa[p_old2]=1/(aa[p_old2]);
      for(m=0;m<n;m++)
          if(m!=l){
 p_old=m*n+l;
 aa[p_old]=aa[p_old]*aa[p_old2];
 }
      for(m=0;m<n;m++){
          for(p=0;p<n;p++){
             if(m!=l){
                 if(p!=l){
```
#define N\_ON N\_dnet+NIZ\*N\_ULD #define NX N\_ON+10 //

\*y=2.0/(1.0+exp(-2.0\*x))-1.0;

a\_r, int\_T a\_c, int\_T b\_c){

for(i=0; i<a\_r ;i++){ for(j=0;j<b\_c;j++){

cc[i\*b\_c+j]=priv;

for(i=0; i<a\_r ;i++){ for(j=0;j<b\_c;j++){

cc[i\*b\_c+j]=priv;

 int\_T i, j, k; real\_T priv;

priv=0;

 int\_T i, j, k; real\_T priv;

priv=0;

 int\_T i, j, k; real\_T priv;

}

 } } }

static

}

 } } }

static

//AxB=C multiplying matrices in line form

for(k=0;k<a\_c;k++) { priv=priv+aa[k+i\*a\_c]\*bb[j+b\_c\*k];

int\_T a\_r, int\_T a\_c, int\_T b\_c, int\_T x\_p){

real\_T \*x = ssGetRealDiscStates(S);

for(k=0;k<a\_c;k++) {

int\_T a\_r, int\_T a\_c, int\_T b\_c, int\_T x\_p){

real\_T \*x = ssGetRealDiscStates(S);

priv=priv+x[x\_p+k+i\*a\_c]\*bb[j+b\_c\*k];

//AxX=C multiplying matrices in line form, x from SimStruc

//XxB=C multiplying matrices in line form, x from SimStruc

#define ETTAQ .051 // Q

#define NDISCSTATES NX

130 MATLAB Applications for the Practical Engineer

static

static

{

}

#define ETTAMI 10.0 //final value of the learning rate

#define ITER\_UCI 20000 //number of steps to stop learning RNN

void ttsig(real\_T x, real\_T \*y) //tansig

void aputab\_c(real\_T \*aa, real\_T \*bb, real\_T \*cc, int\_T

void xputab\_c\_pok(SimStruct \*S, real\_T \*bb, real\_T \*cc,

void aputax\_c\_pok(SimStruct \*S, real\_T \*aa, real\_T \*cc,

```
 pos=m*n+p;
 p_old=m*n+l;
 p_old2=l*n+p;
 aa[pos]=aa[pos]-(aa[p_old]*aa[p_old2]);
 }
 }
 }
 }
 p_old=l*n+l;
      for(p=0;p<n;p++)
          if(p!=l){
 pos=l*n+p;
 aa[pos]=-aa[p_old]*aa[pos];
 }
 }
 pos=n*n;
}
/*====================*
 * S-function methods *
 *====================*/
static
                 void mdlInitializeSizes(SimStruct *S) {
 ssSetNumSFcnParams(S, 0); /* Number of expected parameters */
   if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) {
      return; /* Parameter mismatch will be reported by Simulink */
 }
 ssSetNumContStates(S, 0);
 ssSetNumDiscStates(S, NDISCSTATES);
   if (!ssSetNumInputPorts(S, 1)) return;
 ssSetInputPortWidth(S, 0, NUL+NIZ);
   if (!ssSetNumOutputPorts(S, 2)) return;
 ssSetOutputPortWidth(S, 0, NIZ);
 ssSetOutputPortWidth(S, 1, NIZ*N_ULD);
 ssSetNumSampleTimes(S, 1);
 ssSetNumRWork(S, 0);
 ssSetNumIWork(S, 0);
 ssSetNumPWork(S, 0);
 ssSetNumModes(S, 0);
 ssSetNumNonsampledZCs(S, 0);
 ssSetOptions(S, SS_OPTION_EXCEPTION_FREE_CODE);
}
static
                 void mdlInitializeSampleTimes(SimStruct *S) {
 ssSetSampleTime(S, 0, TS);
 ssSetOffsetTime(S, 0, 0.0);
}
#define MDL_INITIALIZE_CONDITIONS
// Function: mdlInitializeConditions ========================================
static
                 void mdlInitializeConditions(SimStruct *S) {
 real_T *x0 = ssGetRealDiscStates(S);
```
int\_T i, j;

0.6359, 0.32};

 } } } }

static

}

static

int\_T i;

for (i=0;i<NIZ;i++) Y[i]=x[N\_OUT+i];

DY[i]= x[N\_dnet+i];

int\_T i, j, k, priv\_i;

 real\_T WW1[W1][W1\_j]; real\_T WW2[NIZ][W2\_j]; real\_T a\_pr[W1]={0}; real\_T a[W1+1]={0}; real\_T JJ[NIZ]={0}; real\_T W1R[W1\*W1]={0}; real\_T oout[NIZ];

 real\_T WW2pr[NIZ\*W1]={0}; real\_T dadWR[W1\*W1\_ij];

 real\_T dYdW[NIZ\*SIZE\_P]={0}; real\_T dYdW\_T[NIZ\*(SIZE\_P)]; real\_T h\_WxP [NIZ\*SIZE\_P];

 real\_T I\_a[W1]; real\_T PP[W1\*W1\_ij];

#define MDL\_UPDATE

real\_T priv;

for (i=0; i<NIZ\*N\_ULD; i++)

for (i=0; i<W1\_ij+W2\_ij; i++){ x0[i]=.1\*tempRAND[i-(i/20)\*20];} for (i=0; i<SIZE\_P; i++){

> for (j=0; j<SIZE\_P; j++){ if (i==j){

x0[N\_P + i\*SIZE\_P+j]=1000;

 real\_T tempRAND[22]={0.66, 0.01, 0.42, -0.14, -0.39, -0.62, -0.61, 0.36, -0.39, 0.08, -0.69, 0.39, -0.24, 0.72, 0.7, 0.18, -0.01, 0.79, 0.64, 0.28,

Dual Heuristic Neural Programming Controller for Synchronous Generator

http://dx.doi.org/10.5772/58377

133

// Function: mdlOutputs =======================================================

void mdlOutputs(SimStruct \*S, int\_T tid) {

 real\_T \*Y = ssGetOutputPortRealSignal(S, 0); real\_T \*DY = ssGetOutputPortRealSignal(S, 1);

//InputRealPtrsType uPtrs = ssGetInputPortRealSignalPtrs(S,0);

//Function: mdlUpdate ======================================================

void mdlUpdate(SimStruct \*S, int\_T tid) {

real\_T \*x = ssGetRealDiscStates(S);

UNUSED\_ARG(tid); /\* not used in single tasking mode \*/

```
 int_T i, j;
 real_T tempRAND[22]={0.66, 0.01, 0.42, -0.14, -0.39, -0.62, -0.61, 0.36,
0.6359, 0.32};
   for (i=0; i<W1_ij+W2_ij; i++){
 x0[i]=.1*tempRAND[i-(i/20)*20];}
   for (i=0; i<SIZE_P; i++){
       for (j=0; j<SIZE_P; j++){
          if (i==j){
 x0[N_P + i*SIZE_P+j]=1000;
 }
 }
 }
}
// Function: mdlOutputs =======================================================
static
                 void mdlOutputs(SimStruct *S, int_T tid) {
 real_T *Y = ssGetOutputPortRealSignal(S, 0);
 real_T *DY = ssGetOutputPortRealSignal(S, 1);
 real_T *x = ssGetRealDiscStates(S);
 int_T i;
   //InputRealPtrsType uPtrs = ssGetInputPortRealSignalPtrs(S,0);
 UNUSED_ARG(tid); /* not used in single tasking mode */
   for (i=0;i<NIZ;i++)
 Y[i]=x[N_OUT+i];
   for (i=0; i<NIZ*N_ULD; i++)
 DY[i]= x[N_dnet+i];
}
#define MDL_UPDATE
//Function: mdlUpdate ======================================================
static
                 void mdlUpdate(SimStruct *S, int_T tid) {
 int_T i, j, k, priv_i;
 real_T priv;
 real_T WW1[W1][W1_j];
 real_T WW2[NIZ][W2_j];
 real_T a_pr[W1]={0};
 real_T a[W1+1]={0};
 real_T JJ[NIZ]={0};
 real_T W1R[W1*W1]={0};
 real_T oout[NIZ];
 real_T WW2pr[NIZ*W1]={0};
 real_T dadWR[W1*W1_ij];
 real_T I_a[W1];
 real_T PP[W1*W1_ij];
 real_T dYdW[NIZ*SIZE_P]={0};
 real_T dYdW_T[NIZ*(SIZE_P)];
 real_T h_WxP [NIZ*SIZE_P];
```
 pos=m\*n+p; p\_old=m\*n+l; p\_old2=l\*n+p;

> for(p=0;p<n;p++) if(p!=l){

aa[pos]=-aa[p\_old]\*aa[pos];

 } } } }

132 MATLAB Applications for the Practical Engineer

p\_old=l\*n+l;

pos=l\*n+p;

/\*====================\* \* S-function methods \* \*====================\*/

ssSetNumContStates(S, 0);

 ssSetNumSampleTimes(S, 1); ssSetNumRWork(S, 0); ssSetNumIWork(S, 0); ssSetNumPWork(S, 0); ssSetNumModes(S, 0);

ssSetNumNonsampledZCs(S, 0);

 ssSetSampleTime(S, 0, TS); ssSetOffsetTime(S, 0, 0.0);

#define MDL\_INITIALIZE\_CONDITIONS

real\_T \*x0 = ssGetRealDiscStates(S);

 ssSetNumDiscStates(S, NDISCSTATES); if (!ssSetNumInputPorts(S, 1)) return; ssSetInputPortWidth(S, 0, NUL+NIZ);

 ssSetOutputPortWidth(S, 0, NIZ); ssSetOutputPortWidth(S, 1, NIZ\*N\_ULD);

if (!ssSetNumOutputPorts(S, 2)) return;

ssSetOptions(S, SS\_OPTION\_EXCEPTION\_FREE\_CODE);

}

 } pos=n\*n;

static

}

} static

}

static

}

aa[pos]=aa[pos]-(aa[p\_old]\*aa[p\_old2]);

void mdlInitializeSizes(SimStruct \*S) {

return; /\* Parameter mismatch will be reported by Simulink \*/

void mdlInitializeSampleTimes(SimStruct \*S) {

void mdlInitializeConditions(SimStruct \*S) {

// Function: mdlInitializeConditions ========================================

 ssSetNumSFcnParams(S, 0); /\* Number of expected parameters \*/ if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) {

```
 real_T K[SIZE_P*NIZ];
 real_T Ppr[SIZE_P*SIZE_P];
 real_T Kxh_W[SIZE_P*SIZE_P];
 real_T pu[W1_j];
 real_T Apr[NIZ*NIZ]={0};
 real_T Pxh_W[SIZE_P*NIZ] ;
 real_T NNx[SIZE_P] ;
 real_T dnet[NIZ*N_ULD]={0};
 real_T etttami=0;
 int_T SIZE_PxSIZE_P=SIZE_P*SIZE_P;
 int_T W1xW1_ij=W1*W1_ij;
 real_T *x = ssGetRealDiscStates(S);
 InputRealPtrsType uPtrs = ssGetInputPortRealSignalPtrs(S, 0);
 x[N_ON]=x[N_ON]+TS;
 pu[NUL+W1]=1.0;
   for (i=0; i<NUL; i++){
 pu[i]=*uPtrs[i];
 }
   for (i=0; i<W1; i++){
 pu[NUL+i]=x[NH_1+i];
 }
 pu[NUL+W1]=1.0;
   for (i=0; i<NIZ; i++){
 JJ[i]=*uPtrs[NUL+i];
 }
   for (i=0; i<W1; i++) {
      for (j=0; j<W1_j; j++) {
 WW1[i][j]=x[N_W1+i*W1_j+j];
 }
 }
   for (i=0; i<NIZ; i++) {
      for (j=0; j<W2_j; j++) {
 WW2[i][j]=x[N_W2+i*W2_j+j];
 }
 }
   for (i=0; i<W1; i++) {
      for (j=0; j<W1; j++) {
 W1R[i*W1+j]=WW1[i][NUL+j];
 }
 }
 xputab_c_pok(S, pu, a_pr, W1, W1_j, 1, 0);//xputab_c_pok(S,bb, cc, a_r,
a_c, b_c, x_p)
   for (i=0; i<W1; i++) { //a
 ttsig(a_pr[i], &a[i]); // ttsig(real_T x,real_T *y)
 x[NH_1+i]=a[i];
 }
 a[W1]=1.0;
   for (i=0; i<W1; i++)
 I_a[i]=1.0-a[i]*a[i];
 xputab_c_pok(S, a, oout, NIZ, W2_j, 1, N_W2); //Y=WW2*[a;1];
   for (i=0; i<NIZ; i++){
 x[N_OUT+i]=oout[i];
```
JJ[i]=JJ[i]-oout[i];

for (i=0; i<W1\_j; i++)

for (i=1; i<W1; i++){

for (i=0; i<NIZ; i++){

dnet[i\*N\_ULD+j]=0.0;

for (i=0; i<W1\*W1\_ij; i++)

for (j=0; j<W1\_j; j++){ PP[i\*(W1\_ij+W1\_j)+j]=pu[j];

for (j=0; j<N\_ULD; j++){

x[N\_dnet+i\*N\_ULD+j]=priv;

for (i=0; i<W1xW1\_ij; i++) PP[i]=PP[i]+0.9\*dadWR[i]; for(i=0;i<W1;i++){

for(i=0;i<NIZ;i++){

priv=0.0;

for(j=0;j<W1\_ij;j++){

for(j=0;j<W1\_ij;j++){

dYdW[i\*SIZE\_P+j]=priv;

for(j=0; j<W2\_j; j++){

if ( x[N\_ON] > ITER\_UCI ) etttami=1e90\*ETTAMI;

for (i=0; i<NIZ\*NIZ; i=i+NIZ+1){

dYdW[W1\_ij+i\*(W1\_ij+W2\_ij+W2\_j)+j]=a[j];

 transp(dYdW, dYdW\_T, NIZ, SIZE\_P); //dYdW\_T=dYdW' aputax\_c\_pok(S, dYdW, h\_WxP, NIZ, SIZE\_P, SIZE\_P, N\_P);

aputab\_c(h\_WxP, dYdW\_T, Apr, NIZ, SIZE\_P, NIZ);

else etttami=ETTAMI+1000\*ETTAMI/(x[N\_ON]+1.0);

for(k=0;k<W1;k++){

for (k=0; k<W1; k++){

aputax\_c\_pok(S, W1R, dadWR, W1, W1, W1\_ij, N\_dadWR);

x[N\_dadWR+i\*W1\_ij+j]=I\_a[i]\*PP[i\*W1\_ij+j];

priv=priv+x[N\_W2+i\*W2\_j+k]\*x[N\_dadWR+k\*W1\_ij+j];

for(i=0; i<NIZ; i++) { // dydW2

priv=priv+WW2[i][k]\*(1.0-(a[k]\*a[k]))\*WW1[k][j\*N\_TD];}

Dual Heuristic Neural Programming Controller for Synchronous Generator

http://dx.doi.org/10.5772/58377

135

}

 } }

 } }

 } }

}

 } }

 } }

//Kalman

//dnet

PP[i]=0;

PP[i]=pu[i];

priv=0.0;

Dual Heuristic Neural Programming Controller for Synchronous Generator http://dx.doi.org/10.5772/58377 135

```
 JJ[i]=JJ[i]-oout[i];
 }
   for (i=0; i<W1*W1_ij; i++)
 PP[i]=0;
   for (i=0; i<W1_j; i++)
 PP[i]=pu[i];
   for (i=1; i<W1; i++){
      for (j=0; j<W1_j; j++){
 PP[i*(W1_ij+W1_j)+j]=pu[j];
 }
 }
   //dnet
   for (i=0; i<NIZ; i++){
      for (j=0; j<N_ULD; j++){
 priv=0.0;
 dnet[i*N_ULD+j]=0.0;
          for (k=0; k<W1; k++){
 priv=priv+WW2[i][k]*(1.0-(a[k]*a[k]))*WW1[k][j*N_TD];}
 x[N_dnet+i*N_ULD+j]=priv;
 }
 }
 aputax_c_pok(S, W1R, dadWR, W1, W1, W1_ij, N_dadWR);
   for (i=0; i<W1xW1_ij; i++)
 PP[i]=PP[i]+0.9*dadWR[i];
   for(i=0;i<W1;i++){
      for(j=0;j<W1_ij;j++){
 x[N_dadWR+i*W1_ij+j]=I_a[i]*PP[i*W1_ij+j];
 }
 }
   for(i=0;i<NIZ;i++){
      for(j=0;j<W1_ij;j++){
 priv=0.0;
          for(k=0;k<W1;k++){
 priv=priv+x[N_W2+i*W2_j+k]*x[N_dadWR+k*W1_ij+j];
 }
 dYdW[i*SIZE_P+j]=priv;
 }
 }
   for(i=0; i<NIZ; i++) { // dydW2
      for(j=0; j<W2_j; j++){
 dYdW[W1_ij+i*(W1_ij+W2_ij+W2_j)+j]=a[j];
 }
 }
   //Kalman
 transp(dYdW, dYdW_T, NIZ, SIZE_P); //dYdW_T=dYdW'
 aputax_c_pok(S, dYdW, h_WxP, NIZ, SIZE_P, SIZE_P, N_P);
 aputab_c(h_WxP, dYdW_T, Apr, NIZ, SIZE_P, NIZ);
   if ( x[N_ON] > ITER_UCI )
 etttami=1e90*ETTAMI;
   else etttami=ETTAMI+1000*ETTAMI/(x[N_ON]+1.0);
   for (i=0; i<NIZ*NIZ; i=i+NIZ+1){
```
 real\_T K[SIZE\_P\*NIZ]; real\_T Ppr[SIZE\_P\*SIZE\_P]; real\_T Kxh\_W[SIZE\_P\*SIZE\_P];

 real\_T Apr[NIZ\*NIZ]={0}; real\_T Pxh\_W[SIZE\_P\*NIZ] ; real\_T NNx[SIZE\_P] ; real\_T dnet[NIZ\*N\_ULD]={0};

int\_T W1xW1\_ij=W1\*W1\_ij;

for (i=0; i<NUL; i++){ pu[i]=\*uPtrs[i];

for (i=0; i<W1; i++){ pu[NUL+i]=x[NH\_1+i];

for (i=0; i<NIZ; i++){ JJ[i]=\*uPtrs[NUL+i];

for (i=0; i<W1; i++) {

for (i=0; i<NIZ; i++) {

for (i=0; i<W1; i++) {

for (j=0; j<W1\_j; j++) { WW1[i][j]=x[N\_W1+i\*W1\_j+j];

for (j=0; j<W2\_j; j++) { WW2[i][j]=x[N\_W2+i\*W2\_j+j];

for (j=0; j<W1; j++) { W1R[i\*W1+j]=WW1[i][NUL+j];

for (i=0; i<W1; i++) { //a

int\_T SIZE\_PxSIZE\_P=SIZE\_P\*SIZE\_P;

real\_T \*x = ssGetRealDiscStates(S);

InputRealPtrsType uPtrs = ssGetInputPortRealSignalPtrs(S, 0);

xputab\_c\_pok(S, pu, a\_pr, W1, W1\_j, 1, 0);//xputab\_c\_pok(S,bb, cc, a\_r,

xputab\_c\_pok(S, a, oout, NIZ, W2\_j, 1, N\_W2); //Y=WW2\*[a;1];

ttsig(a\_pr[i], &a[i]); // ttsig(real\_T x,real\_T \*y)

real\_T pu[W1\_j];

134 MATLAB Applications for the Practical Engineer

real\_T etttami=0;

 x[N\_ON]=x[N\_ON]+TS; pu[NUL+W1]=1.0;

pu[NUL+W1]=1.0;

}

}

}

 } }

 } }

 } }

}

a\_c, b\_c, x\_p)

a[W1]=1.0;

x[NH\_1+i]=a[i];

for (i=0; i<W1; i++) I\_a[i]=1.0-a[i]\*a[i];

for (i=0; i<NIZ; i++){ x[N\_OUT+i]=oout[i];

```
 Apr[i]=Apr[i]+etttami;
 }
 invmatrix(Apr, NIZ);
 xputab_c_pok(S, dYdW_T, Pxh_W, SIZE_P, SIZE_P, NIZ, N_P);
 aputab_c(Pxh_W, Apr, K, SIZE_P, NIZ, NIZ); //K=PM*h_WM*A;
 aputab_c(K, JJ, NNx, SIZE_P, NIZ, 1);
   for(i=0; i<SIZE_P; i++){
 x[i]=x[i] +NNx[i];
 }
 aputab_c(K, dYdW, Kxh_W, SIZE_P, NIZ, SIZE_P); //K*h_WM
 aputax_c_pok(S, Kxh_W, Ppr, SIZE_P, SIZE_P, SIZE_P, N_P);
    for (i=0; i<SIZE_P; i++){ Ppr[i*(SIZE_P+1)]=Ppr[i*(SIZE_P+1)] - ETTAQ/
(x[N_ON]+1.0);
 }
   for (i=0; i<SIZE_PxSIZE_P; i++)
 x[N_P+i]=x[N_P+i]-1*Ppr[i];
 UNUSED_ARG(tid);
}
// mdlTerminate
static
                  void mdlTerminate(SimStruct *S) {
 int_T i, j;
 real_T *x = ssGetRealDiscStates(S);
 printf("\n");
   for(i=0; i<W1*W1_j+NIZ*W2_j; i++){
 printf("%f\t", x[i]);
 }
 UNUSED_ARG(S);
}
#ifdef MATLAB_MEX_FILE /* Is this file being compiled as a MEX-file? */
#include "simulink.c" /* MEX-file interface mechanism */
#else
#include "cg_sfun.h" /* Code generation registration function */
#endif
```
**References**

2004.

1994.

5(14), 279-297.

trand Reinhold, pp493-525.

997-1007, Sept 1997.

[1] P.Werbos "Aproximate dynamic programming for real-time control and neural mod‐ eling" in Handbook of Intelligent Control, White and Sofge Eds. New York: Van Nos‐

Dual Heuristic Neural Programming Controller for Synchronous Generator

http://dx.doi.org/10.5772/58377

137

[2] D. Prokhorov: "Adaptive critic design", IEEE Trans. Neural Networks, vol. 8, pp.

[3] Wiliams, R. J. and Zisper D: "A learning agorithm for continually running fully re‐

[4] Pedro DeLima: "Application of Adaptive Critic Designs for Fault Tolerant Control", IEEE Computational intelligence Society Walter J Karplus Summer Research Grant

[5] R.E. Kalman, ''A new approach to linear filtering and prediction problems,''Transac‐

[6] G.V. Puskorius and L. A. Feldkamp: "Neurocontrol of Nonlinear Dynamic Systems with Kalman Filter Trained Recurrent Networks", IEEE Trans. on Neural Networks,

[8] M. Cernasky, L. Benuskova:" Simple recurrent network trained by RTRL and extend‐

[9] Narendra K.S., Parthasarathy K. "Identification and Control of Dynamical Systems Using Neural Networks", IEEE transaction on Neural Networks, 1:4-27, March 1990.

[10] IEEE recommended practice for excitation system models for power system stability

[11] W. Mielczarski and A. M. Zajaczkowski "Nonlinear Field Voltage Control of a Syn‐ chronous Generator using Feedback Linearization", Automatica Vol 30, pp1625-1630,

[12] P. Shamsollahi, O. P. Malik, "Real time Implementation and Experimental Studies of a Neural Adaptive Power System Stabilizer", IEEE Trans. on Energy Conversion, Vol.

[13] G. K. Venayagamoorthy, Ronald G. Harley, and Donald Wuncsh "Implementation off an adaptive neural network for efective control of turbogenerators" *IEEE Budapest*

[14] T. Kobayashi and A. Yokoyama, "An Adaptive Neuro-Control System of Synchro‐ nous Generator for Power system stabilization", IEEE Trasaction on Energy Conver‐

ed Kalman filter algorithms" Neural Network World 13(3) (2003) 223–234

tions of the ASME, Ser. D, Journal of Basic Engineering, 82, 35–45 (1960).

current neural networks", Neural Copmutation, 1:270-280.

[7] Simon S. Haykin: "Kalman Filtering and Neural Networks", 2001

studies, : IEEE St. 421.5-2002 (Section 9).

*Power Tech. Conf. 1999. paper BPT99 pp 431 – 23.*

sion, Vol. 11, No. 3, September 1996.

14. No. 3, September 1999.

#### **Author details**

Mato Miskovic1 , Ivan Miskovic2 and Marija Mirosevic3

1 HEP Hydro-power Dubrovnik, Dubrovnik, Croatia

2 Brodarski Institute, Zagreb, Croatia

3 University of Dubrovnik, Dubrovnik, Croatia

#### **References**

Apr[i]=Apr[i]+etttami;

for(i=0; i<SIZE\_P; i++){ x[i]=x[i] +NNx[i];

aputab\_c(K, JJ, NNx, SIZE\_P, NIZ, 1);

for (i=0; i<SIZE\_PxSIZE\_P; i++) x[N\_P+i]=x[N\_P+i]-1\*Ppr[i];

for(i=0; i<W1\*W1\_j+NIZ\*W2\_j; i++){

 xputab\_c\_pok(S, dYdW\_T, Pxh\_W, SIZE\_P, SIZE\_P, NIZ, N\_P); aputab\_c(Pxh\_W, Apr, K, SIZE\_P, NIZ, NIZ); //K=PM\*h\_WM\*A;

 aputab\_c(K, dYdW, Kxh\_W, SIZE\_P, NIZ, SIZE\_P); //K\*h\_WM aputax\_c\_pok(S, Kxh\_W, Ppr, SIZE\_P, SIZE\_P, SIZE\_P, N\_P);

real\_T \*x = ssGetRealDiscStates(S);

#include "simulink.c" /\* MEX-file interface mechanism \*/

and Marija Mirosevic3

for (i=0; i<SIZE\_P; i++){ Ppr[i\*(SIZE\_P+1)]=Ppr[i\*(SIZE\_P+1)] - ETTAQ/

void mdlTerminate(SimStruct \*S) {

#ifdef MATLAB\_MEX\_FILE /\* Is this file being compiled as a MEX-file? \*/

#include "cg\_sfun.h" /\* Code generation registration function \*/

invmatrix(Apr, NIZ);

}

136 MATLAB Applications for the Practical Engineer

}

}

static

}

#else

#endif

**Author details**

Mato Miskovic1

}

}

(x[N\_ON]+1.0);

// mdlTerminate

int\_T i, j;

printf("\n");

UNUSED\_ARG(S);

printf("%f\t", x[i]);

, Ivan Miskovic2

3 University of Dubrovnik, Dubrovnik, Croatia

2 Brodarski Institute, Zagreb, Croatia

1 HEP Hydro-power Dubrovnik, Dubrovnik, Croatia

UNUSED\_ARG(tid);

[15] D.Flynn, S. McLonne, G. W. Irwin, M. D. Brown, E. Swidenbank, and B. W. Hogg, "Neural control turbogeneraroe systems", Automatica, vol 33, no,11, pp. 1961 – 1973, 1997.

**Chapter 5**

**Design of Fractionation Columns**

Additional information is available at the end of the chapter

Distillation is a physical process used to separate a fluid mixture of two (binary) or more (multicomponent) substances into its component parts. In most cases, the components to be separated are miscible liquids with different volatilities and boiling points. This separation process is a thermal unit operation that utilizes the differences of vapour pressure to produce the separa‐ tion. In this process, the vapour or liquid mixture is heated whereby the more volatile components are evaporated, condensed and allowed to drip or drip apart, i.e. distil or destillare, as it was originally called in Latin. It is from this fact of 'Dripping' that the name of

Distillation may be carried out in an intermittent or batch process or it may be carried out in a continuous steady-state process with a continuous feed stream. Distillation systems may also operate at different pressures, with higher pressures used for the separation of the more volatile materials and reduced pressures, in vacuum distillation systems, for heavier higher

Simple distillation of alcoholic beverages and petroleum has been practiced for a very long time. In the Middle Ages, Damascus played a part as a distilling centre, where distillation came into its own as the best and quickest way of obtaining pure chemical substances [1]. In the writings of the thirteenth century Syrian Scholar, Al-Dimashki, we find the earliest reference to the industrial distillation of crude oil or petroleum. Al-Dimashki says: "Many types of petroleum are water white by nature and so volatile that they can not be stored in open vessels. Others are obtained from a kind of pitch or bitumen in a turbid and dark condition, but by further treatment they can be made clear and white by distilling them like

An earlier account of crude oil distillation was given by Rhases (865-925) in his book Sirr alasrar (Secretum secretorum or Secret of secrets), and in 1512 Hieronymus Brunschwig wrote a book

> © 2014 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Hassan Al-Haj Ibrahim

http://dx.doi.org/10.5772/57249

the distillation process was derived.

boiling-point materials, lowering thereby their boiling points.

**1. Introduction**

rose-water" [1].

**Chapter 5**