**Tool of the Complete Optimal Control for Variable Speed Electrical Drives**

Marian Gaiceanu

Proceedings of the International Conference of Aerospace Sciences "AEROSPATIAL

[17] Grigorie, T. L., Botez, R. M., Lungu, M., Edu, I. R. & Obreja, R. (2012 a). MEMS Gyro Performance Improvement through Bias Correction over Temperature using an Adaptive Neural Network trained Fuzzy Inference System, *Proceedings of the Institu‐ tion of Mechanical Engineers, Part G: Journal of Aerospace Engineering*, September 2012,

[18] Grigorie, T. L., Obreja, R., Sandu, D. G., Corcau, J. I. (2012 b). Allan variance analysis of the miniaturized sensors in a strap-down inertial measurement unit, 12th Interna‐ tional Multidisciplinary Scientific GeoConference - SGEM2012, ISSN 1314-2704, Vol.

[19] Hopkins, R., Barbour, N., Gustafson, D.E. & Sherman, P. (2010). Miniature Inertial and Augmentation Sensors for Integrated Inertial/GPS Based Navigation Applica‐ tions, NATO RTO Lecture Series, RTO-EN-SET-116, Low-Cost Navigation Sensors

[20] Kraft, M. (2000). Micromachined inertial sensors: The state-of-the-art and a look into

[21] KVH Industries Inc. (2007). *An update on KVH fiber optic gyros and their benefits relative*

[22] Lawrence, A. (1998). *Modern Inertial Technology: Navigation, Guidance and Control* – Second Edition, Springer Verlag, New York, ISBN: 978-1-4612-7258-8, 1998

[23] Pavlath, G. (2006). Fiber Optic Gyros: The Vision Realized, 18th International Confer‐

[24] Radix, J.C. (1993). *Systemes inertiels a composants lies <<Strap-Down>>, Cepadues-Edi‐ tions*, Ecole Nationale Superieure de l'Aeronautique et de l'Espace SUP'AERO, Tou‐

[25] Salychev, O. S. (1998). *Inertial Systems in Navigation and Geophysics*, Bauman MSTU

[26] Savage, P. G. (2000). *Strapdown Analytics, Part 1.* Strapdown Associates, Inc., ISBN:

[27] Schmidt, G. (2010). INS/GPS Technology Trends, NATO RTO Lecture Series, RTO-EN-SET-116, Low-Cost Navigation Sensors and Integration Technology, March 2010

[28] Tawney, J., et al. (2006). Photonic Crystal Fiber IFOGs, Optical Society of America 18th International Conference on Optical Fiber Sensors, Cancun, Mexico, October

[29] Titterton, D. H. & Weston, J. (2004). *Strapdown inertial navigation technology - 2nd Edi‐ tion*, Institution of Engineering and Technology, ISBN: 0 863413587, USA, 2004.

2010", Section 5, pp. 1-10, ISSN 2067-8622, Bucharest, 20-21 October, 2010

Vol. 226(9), pp. 1121-1138, ISSN: 0954-4100, 2012

3, pp. 443 – 450, June 17-23, 2012

338 MATLAB Applications for the Practical Engineer

and Integration Technology, March 2010

*to other gyro technologies*, March 2007

louse, France

2006

9780971778603, USA, 2000

the future, *Meas. Control*, vol. 33, no. 6, pp.164 -168, 2000

ence on Optical Fiber Sensors, Cancun, Mexico, October 2006

Press, ISBN: 5703813468, Moscow, Russia, 1998

Additional information is available at the end of the chapter

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

#### **1. Introduction**

The main objective of the chapter is to use Matlab software in order to develop the optimal control knowledge in relation to the biggest worldwide power consumers: electrical machines. These two major components, optimal control and electrical machines form an optimal drive system. The proposed drive system conducts to the energy efficiency increasing during transient regimes (starting, stopping and reversing), therefore reducing the impact upon the environment. Taking into account the world energy policies, the chapter is strategically oriented towards the compatibility with the priority requirements from world research programmes. This chapter will offer an original Matlab tool in order to implement a highly performant control for electrical drives. The optimal solution provided by the optimal problem is obtained by numerical integration of the matrix Riccati differential equation (MRDE). The proposed optimal control, based on energetic criterion, can be implemented on-line by using a real time experimental platform. The solution has three terms: the first term includes the reference state of the electrical drive system, the second term assures a fast compensation of the disturbance (i.e. the load torque in case of the electrical drives) by using feedforward control, and the third is the state feedback component. Therefore, there is a completed optimal control for linear dynamic systems. The simulation and experimental results will be provided. By using knowledge regarding the electrical machines, electrical and mechanical measure‐ ments, control theory, digital control, and real time implementation, a high degree of inter‐ disciplinarity is obtained in the developed Matlab tool The optimal problem statement is without constraints; by choosing adequate weighting matrices, the magnitude constraints of control and of state variable are solved.

© 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.

Matlab is a widespread software used both in academia and research centres. The author of this chapter facilitates the easy understanding and implementation of the highly efficient electrical drive systems by using Matlab software. The energy problem is one of the most important aspects related to the sustainable growth of the industrial nations. On the one hand, the industrial systems and the increased quality of life require every year a higher quantity of electrical energy produced; on the other hand, the same quality of life requires the generation and the harnessing of this energy with low pollution.

optimization based on the speed controller which can be optimized in the time or in the frequency range. In the time range, the symmetrical optimum method can be used. In the frequency range, the optimizing is done according to the amplitude and phase margin rules, known from the control theory, but these are the conventional controls that are used in drives.

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

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

341

The interdisciplinarity of the chapter consists of using specific knowledge from the fields of: energy conversion, power converters, Matlab/Simulink simulation software, real time implementation based on dSPACE platform, electrotechnics, and advanced control techni‐

This chapter is focused on mathematical modelling, optimal control development on power inverter, Matlab implementation and analysis of the optimal control solution in order to minimize the electric input energy of the drive system. The Matlab on-line solution can be downloaded in real-time platforms, as dSPACE or dsPIC Digital Signal Controllers (Gaiceanu,

The practical implication of using the proposed method is the on-line calculation of the optimal

The strategy of the optimal control electrical drive system is to use both the conventional control during the steady state and the developed optimal control during the transients. The conventional control is based on the rotor field orientated control using Proportional Integral controllers. If the speed error is different from zero, the software switch will enable the optimal control, else will enable the conventional control. Both types of control are developed in the

**2. Mathematical model of the variable speed electrical drives (VSED)**

There are variable speed DC and AC drives. A DC drive system controlled by armature voltage at constant field is a linear, invariant and controllable dynamic system. Variable speed AC drives are mostly used in connection with induction and synchronous electrical machines. The modern drives contain the energy source, the power converters, the electrical machines, the mechanical power transmission and the load, linked by an adequate control. The well-known mathematical model of the vector controlled AC drive using three phase squirrel cage induction machine in rotor field coordinates supplied from the current inverter, operating at the constant flux (Leonhard, 1996) becomes a linear system (Gaiceanu, 2002). For the electrical drives having the power less than 50 kW, the power loss of the converter could be neglected. In order to compute the optimal solution for synchronous Drive system, the control in rotor

control solution. The main lines of the research methodology are:

**•** The optimal control of the dynamical electrical drives problem formulation;

ques.

*et al.*, 2008).

next Sections.

**•** The optimal control solution;

**•** The numerical simulation results;

**•** The implementation of the optimal control.

On the one hand, the electrical machines are the main worldwide energy consumer (Siemens), (Tamimi, *et al.*, 2006). Therefore, the focus of the research on this system type is essential in order to minimize the energetical losses. A seemingly insignificant reduction of electricity consumed by each electrical motor has particularly important consequences at national/ international level, primarily by reducing consumption of fossil fuels widely used to generate electricity, therefore contributing also to CO2 emissions reduction. Researches in this area are at different levels: I) direct conversion by redesigning the electric motors to obtain high efficiency; II) the power electronic systems by introducing static converters with diminished switching losses; III) control implementation of the conversion process by using the enhanced capabilities of new generation digital signal processor DSP (Veerachary M., 2002).

The efficiency of the motor is high on steady state operation, but very low in transients. On the other hand, the associated equipment for energy conversion purposes, i.e. power convert‐ ers together with the control subsystem, assures both high static and dynamic performances (related to the settling time, error minimizing and so on), but they are not related to energy expenditure. The analysis of various technical publications shows a yearly increase of energy expenditure of the electrical motors, with a rate of 1,5% in Europe, 2,2% in USA, 3,1% in Japan; these values can even reach the rate of 13%, depending on the development of calculus technology (Matinfar, *et al.,* 2005), (Gattami, *et al.,*2005), (COM, 2006).

Dynamic programming tool is very useful to obtain the optimal solution, but the values of the torque and speed must be known in advance (Mendes, *et al.*, 1996). In machine building industry the operating cycle is well-known (Lorenz, *et al.*, 1992; 1998) therefore, the dynamic programming is an efficient tool only for these type of applications. But the optimization method fails when the conventional electrical drives are taken into account because the load torque could not be known in advance, i.e. at the final time of the process.

The optimal solution (Gaiceanu, 1997) provided in this chapter overcomes this type of scientific barrier (Rosu, 1985), making possible implementation of the optimal solution to any type of application, without knowing the speed and load torque in advance. As it is nonrecursive, determined on the basis of variational methods (by integrating the Matrix Riccati Differential Equation), the optimal solution can be calculated on-line without memorizing it from the final time to the initial time, as in the recursive method. Therefore, the nonrecursive solution can be implemented to any type of electrical motor with any load variation, for any operating dynamical regimes.

On the other hand, the problem of increasing power conversion efficiency in the drive systems is the subject of numerous worldwide studies. The Siemens Automation (Siemens) claims the optimization based on the speed controller which can be optimized in the time or in the frequency range. In the time range, the symmetrical optimum method can be used. In the frequency range, the optimizing is done according to the amplitude and phase margin rules, known from the control theory, but these are the conventional controls that are used in drives.

The interdisciplinarity of the chapter consists of using specific knowledge from the fields of: energy conversion, power converters, Matlab/Simulink simulation software, real time implementation based on dSPACE platform, electrotechnics, and advanced control techni‐ ques.

This chapter is focused on mathematical modelling, optimal control development on power inverter, Matlab implementation and analysis of the optimal control solution in order to minimize the electric input energy of the drive system. The Matlab on-line solution can be downloaded in real-time platforms, as dSPACE or dsPIC Digital Signal Controllers (Gaiceanu, *et al.*, 2008).

The practical implication of using the proposed method is the on-line calculation of the optimal control solution. The main lines of the research methodology are:


Matlab is a widespread software used both in academia and research centres. The author of this chapter facilitates the easy understanding and implementation of the highly efficient electrical drive systems by using Matlab software. The energy problem is one of the most important aspects related to the sustainable growth of the industrial nations. On the one hand, the industrial systems and the increased quality of life require every year a higher quantity of electrical energy produced; on the other hand, the same quality of life requires the generation

On the one hand, the electrical machines are the main worldwide energy consumer (Siemens), (Tamimi, *et al.*, 2006). Therefore, the focus of the research on this system type is essential in order to minimize the energetical losses. A seemingly insignificant reduction of electricity consumed by each electrical motor has particularly important consequences at national/ international level, primarily by reducing consumption of fossil fuels widely used to generate electricity, therefore contributing also to CO2 emissions reduction. Researches in this area are at different levels: I) direct conversion by redesigning the electric motors to obtain high efficiency; II) the power electronic systems by introducing static converters with diminished switching losses; III) control implementation of the conversion process by using the enhanced

capabilities of new generation digital signal processor DSP (Veerachary M., 2002).

technology (Matinfar, *et al.,* 2005), (Gattami, *et al.,*2005), (COM, 2006).

torque could not be known in advance, i.e. at the final time of the process.

dynamical regimes.

The efficiency of the motor is high on steady state operation, but very low in transients. On the other hand, the associated equipment for energy conversion purposes, i.e. power convert‐ ers together with the control subsystem, assures both high static and dynamic performances (related to the settling time, error minimizing and so on), but they are not related to energy expenditure. The analysis of various technical publications shows a yearly increase of energy expenditure of the electrical motors, with a rate of 1,5% in Europe, 2,2% in USA, 3,1% in Japan; these values can even reach the rate of 13%, depending on the development of calculus

Dynamic programming tool is very useful to obtain the optimal solution, but the values of the torque and speed must be known in advance (Mendes, *et al.*, 1996). In machine building industry the operating cycle is well-known (Lorenz, *et al.*, 1992; 1998) therefore, the dynamic programming is an efficient tool only for these type of applications. But the optimization method fails when the conventional electrical drives are taken into account because the load

The optimal solution (Gaiceanu, 1997) provided in this chapter overcomes this type of scientific barrier (Rosu, 1985), making possible implementation of the optimal solution to any type of application, without knowing the speed and load torque in advance. As it is nonrecursive, determined on the basis of variational methods (by integrating the Matrix Riccati Differential Equation), the optimal solution can be calculated on-line without memorizing it from the final time to the initial time, as in the recursive method. Therefore, the nonrecursive solution can be implemented to any type of electrical motor with any load variation, for any operating

On the other hand, the problem of increasing power conversion efficiency in the drive systems is the subject of numerous worldwide studies. The Siemens Automation (Siemens) claims the

and the harnessing of this energy with low pollution.

340 MATLAB Applications for the Practical Engineer


The strategy of the optimal control electrical drive system is to use both the conventional control during the steady state and the developed optimal control during the transients. The conventional control is based on the rotor field orientated control using Proportional Integral controllers. If the speed error is different from zero, the software switch will enable the optimal control, else will enable the conventional control. Both types of control are developed in the next Sections.

### **2. Mathematical model of the variable speed electrical drives (VSED)**

There are variable speed DC and AC drives. A DC drive system controlled by armature voltage at constant field is a linear, invariant and controllable dynamic system. Variable speed AC drives are mostly used in connection with induction and synchronous electrical machines. The modern drives contain the energy source, the power converters, the electrical machines, the mechanical power transmission and the load, linked by an adequate control. The well-known mathematical model of the vector controlled AC drive using three phase squirrel cage induction machine in rotor field coordinates supplied from the current inverter, operating at the constant flux (Leonhard, 1996) becomes a linear system (Gaiceanu, 2002). For the electrical drives having the power less than 50 kW, the power loss of the converter could be neglected. In order to compute the optimal solution for synchronous Drive system, the control in rotor field oriented reference frame of surface mounted permanent magnet synchronous motor (PMSM) drives is used. The decoupling of the control loops (torque and flux) is performed in rotor field reference frame. By setting a zero value of direct current stator component, the mathematical model of the PMSM becomes linear. Therefore, there is a common mathematical model of the electrical drives, described by the linear state space standard form representation: control loops (torque and flux) is performed in rotor field reference frame. By setting a zero value of direct current stator component the mathematical model of the PMSM becomes linear. Therefore, there is a common mathematical model of the electrical drives, described by the linear state space standard form representation:

$$\dot{\mathbf{x}}(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}u(t) + \mathbf{G}w(t) \tag{1}$$

squirrel cage induction machine in rotor field coordinates supplied from the current inverter, operating at the constant flux (Leonhard, 1996) becomes a linear system (Gaiceanu, 2002). For the electrical drives having the power less than 50 kW, the power loss of the

permanent magnet synchronous motor (PMSM) drives is used. The decoupling of the

e

w

*<sup>i</sup> T ii t i q t p Ti J F T T pt p T ki i*

e e

*r s*

*mm q*

<sup>ï</sup> =- - <sup>ï</sup> <sup>ï</sup> =×× ïî

,

w

– electromagnetic torque

<sup>2</sup> (c) 3 1 *m*

s= × <sup>+</sup>

*r L*

d (a

*qs r mr*

where rotor magnetizing current;

, stator and rotor time constants;

*q w*

mathematical model can be obtained (Fig.2):

(b)

*L*

)

T

e

s Tr

−

e

d

d

d

 

 

, *mr ds*

e

m

(1)


*i ii f*


*q*

with magnetizi g

– electromagnetic

m mr q

=− −

J F T T

e

ω

pt p T ki i

T ii t i q t p Ti

e

=⋅⋅

ω

= +

d d

i

mr r mr ds

+ =

*s*

angle of the magnetizing current vector; , electrical and mechanical rotor speed

etizing inductance, rotor magnetizing dispersion factor.

mL

<sup>2</sup> (c) 3 1

r

*r*


<sup>σ</sup> = ⋅ <sup>+</sup>

% INITIAL DATA isd=Idref; imr=isd;

Fv=0.0006;

 G=[-1/J; 0 ]; % The output matrix Ct=[1 0]; % D=[0 0];

**Figure 2.** Standard space state form characterization of the three phase Induction Motor

m

−

L

p

ω ω

=

*s*

K p

*e*

*w*

By varying the flux component of the stator current, *ids*, the magnetizing current response, *im*r(*t*), can have signifiant delay. By using an adequate control of power inverter (Gaiceanu, 2002), the rotor magnetizing current, *i*mr, can be maintained at the constant value, and according to the mathematical model of the IM (1a), the flux component of the stator current, *ids*, equates the rotor magnetizing current): *i*mr=*ids*=ct.. Therefore, the standard state space form of the IM

km=2\*Lm\*imr/(3\*(1+sigmar));

% The matrix of the system A=[-Fv/J 0; 1 0]; % The control vector kq=1/(Tr\*imr); B=[ km/J; kq ];

% The perturbation vector

*TL p F*

(


T T

e e

,

ω

s

mr ds

i ii f

, stator and rotor time constants;

b)

n

− −

d (a)

qs r mr

, lux and torque components of the stator current;

L

,

angle of the magnetizing current vector; , electrical and mechanical rotor speed,

e

−

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

w

inductance, rotor magnetizing dispersion factor

By varying the flux component of the stator current, ids, the magnetizing current response, im<sup>r</sup>(t), can have signifiant delay. By using an adequate control of power inverter (Gaiceanu, 2002), the rotor magnetizing current, imr, can be maintained at the constant value, and according to the mathematical model of the IM (1a), the flux component of the stator current, ids, equates the rotor magnetizing current): imr=ids=ct.. Therefore, the standard state

−

rotor magnetizing current; , lux and torque components of the

o

s

q

r

−

s

space form of the IM mathematical model can be obtained (Fig.2):

Figure 2. Standard space state form characterization of the three phase Induction Motor.

– l

(2)

ad torque, number of pole pairs, – viscous force;

p F

343

stator current;

.

– load torque, number of pole pairs, – viscous force;

q w

torque,

L

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

e

w

d d d

ì

ï ï ï = + ï í ï

*mr r mr ds*

ï + =

e

*e*

*Ts Tr*

e

*m*

*K p*

w w

with magn


*m*

*L*

*p*

=

*T*

d d

where *x*(*t*) is the state vector, **u**(t) is the control vector, and the perturbation vector is the load torque w(*t*)= *TL* (*t*) . where x t)( is the state vector, u(t) is the control vector and the perturbation vector is the load torque w(t)= [ tT )( ] <sup>L</sup> .

Taking into consideration the most used electrical machine, the three-phase squirrel cage induction machine, the electrical drive based on it will be explained as follows. Starting from the nameplate data of the induction machine (Fig.1), the main parameters of the electrical motor can be found (Appendix 1). Taking into consideration the most used electrical machine, three-phase squirrel cage induction machine, the electrical drive based on it will be explained as follows. Starting from the nameplate data of the induction machine (Fig.1), the main parameters of the

electrical motor can be found (Appendix 1).

**Figure 1.** Nameplate data of the electrical machine (three phase induction motor)

Taking into consideration the mathematical model of the three phase squirrel cage induction machine (IM) in rotor field reference system operating at constant flux (Leonhard, 1996): The mathematical model of the three phase squirrel cage induction machine (IM) in rotor field reference system operating at constant flux is (Leonhard, 1996):

Figure 1. Nameplate data of the electrical machine (three phase induction motor)

torque,

q w

rotor magnetizing current; , lux and torque components of the

angle of the magnetizing current vector; , electrical and mechanical rotor speed,

e

−

w

inductance, rotor magnetizing dispersion factor

−

o

s

q

r

−

s

– l

L

stator current;

.

ad torque, number of pole pairs, – viscous force;

p F

$$\begin{cases} \Gamma\_r \frac{d\dot{l}\_{mr}}{dt} + \dot{l}\_{mr} = \dot{l}\_{dc} \\ \frac{d\dot{q}\_1}{dt} = \frac{q\_2}{p} + \frac{\dot{l}\_{qr}}{\Gamma\_r \dot{l}\_{mr}} \\ \frac{d}{p} \frac{d\dot{q}\_2}{dt} = \Gamma\_r - \frac{\dot{F}}{p} \mathbf{q}\_1 - \Gamma\_L \\ \Gamma\_s = k\_{sr} \cdot \dot{l}\_{mr} - \dot{l}\_{pr} \end{cases} (a)$$

field oriented reference frame of surface mounted permanent magnet synchronous motor (PMSM) drives is used. The decoupling of the control loops (torque and flux) is performed in rotor field reference frame. By setting a zero value of direct current stator component, the mathematical model of the PMSM becomes linear. Therefore, there is a common mathematical model of the electrical drives, described by the linear state space standard form representation:

where *x*(*t*) is the state vector, **u**(t) is the control vector, and the perturbation vector is the load

load torque w(t)= [ tT )( ] <sup>L</sup> .

Taking into consideration the most used electrical machine, the three-phase squirrel cage induction machine, the electrical drive based on it will be explained as follows. Starting from the nameplate data of the induction machine (Fig.1), the main parameters of the electrical

%Nominal Power (mechanical power)

%total equivalent inertia mass J=0.002; %[kgm^2] % number of pole pairs

The mathematical model of the three phase squirrel cage induction machine (IM) in rotor field

 Pn=750;% kW %Stator resistance Rs=1.7; %[ohmi] %rotor resistance Rr=2.55; %[ohmi] %stator leakage inductance Lss=0.00986; %[H] %rotor leakage inductance Lrs=0.01002; %[H] %mutual inductance Lm=0.2680; %[H]

 p=2; %power factor cosfi=0.71; %rated efficiency etan=0.794; %rated speed nn=1480; %frequency f1=50; %[Hz]

**Figure 1.** Nameplate data of the electrical machine (three phase induction motor)

reference system operating at constant flux is (Leonhard, 1996):

electrical motor can be found (Appendix 1).

torque w(*t*)= *TL* (*t*) .

motor can be found (Appendix 1).

342 MATLAB Applications for the Practical Engineer

*x Ax Bu Gw* &() () () () *t tt t* = ++ (1)

Figure 1. Nameplate data of the electrical machine (three phase induction motor)

Taking into consideration the mathematical model of the three phase squirrel cage induction

by the linear state space standard form representation:

xɺ t)( = Ax t)( + Bu t)( + Gw t)( (1)

squirrel cage induction machine in rotor field coordinates supplied from the current inverter, operating at the constant flux (Leonhard, 1996) becomes a linear system (Gaiceanu, 2002). For the electrical drives having the power less than 50 kW, the power loss of the converter could be neglected. In order to compute the optimal solution for synchronous Drive system, the control in rotor field oriented reference frame of surface mounted permanent magnet synchronous motor (PMSM) drives is used. The decoupling of the

where x t)( is the state vector, u(t) is the control vector and the perturbation vector is the where rotor magnetizing current; , lux and torque components of the stator current; – electromagnetic torque , *mr ds q e s i ii f T* - - – load torque, number of pole pairs, – viscous force; , stator and rotor time constants; *TL p F Ts Tr* - (2) e b) ( p ω ω =

n

 angle of the magnetizing current vector; , electrical and mechanical rotor speed , *e w q w* - -

r

with magnetizi g

– electromagnetic

m mr q

$$
\alpha \rho\_e = p \alpha \rho \tag{5}
$$

from the nameplate data of the induction machine (Fig.1), the main parameters of the with magn *m L* etizing inductance, rotor magnetizing dispersion factor. *r s* - <sup>2</sup> (c) 3 1 mL <sup>σ</sup> = ⋅ <sup>+</sup>

m

−

L

e

pt p T ki i

d (a)

qs r mr

,

s

mr ds

− −

, stator and rotor time constants;

d d

i

= +

mr r mr ds

+ =

t i q t p Ti

ω

e

e

s Tr

−

m

(1)

K p

d

d

 

$$K\_m = \frac{2}{3}p \cdot \frac{L\_m}{1 + \sigma\_r} \qquad \text{ (c)}$$

By varying the flux component of the stator current, *ids*, the magnetizing current response, *im*r(*t*), can have signifiant delay. By using an adequate control of power inverter (Gaiceanu, 2002), the rotor magnetizing current, *i*mr, can be maintained at the constant value, and according to the mathematical model of the IM (1a), the flux component of the stator current, *ids*, equates the rotor magnetizing current): *i*mr=*ids*=ct.. Therefore, the standard state space form of the IM mathematical model can be obtained (Fig.2): By varying the flux component of the stator current, ids, the magnetizing current response, im<sup>r</sup>(t), can have signifiant delay. By using an adequate control of power inverter (Gaiceanu, 2002), the rotor magnetizing current, imr, can be maintained at the constant value, and according to the mathematical model of the IM (1a), the flux component of the stator current, ids, equates the rotor magnetizing current): imr=ids=ct.. Therefore, the standard state space form of the IM mathematical model can be obtained (Fig.2):

$$\begin{array}{l} \v{} \quad \text{TMTTA} \quad \text{DATA} \\ \quad \text{ids-I-derf;} \\ \quad \text{imr-1-id} \\ \quad \text{km-2} \,\text{Lm-1-in} \\ \quad \text{Fv-0.000;} \\ \quad \text{\(\uparrow\)} \\ \quad \text{The matrix of the system} \\ \quad \text{\(\uparrow\)} \\ \quad \text{\(\uparrow\)} \\ \quad \text{The control vector} \\ \quad \text{\(\uparrow\)} \\ \quad \text{\(\uparrow\)} \\ \quad \text{\(\uparrow\)} \\ \quad \text{\(\uparrow\)} \\ \quad \text{\(\uparrow\)} \\ \quad \text{\(\uparrow\)} \\ \quad \text{\(\uparrow\)} \\ \quad \text{\(\uparrow\)} \\ \quad \text{\(\uparrow\)} \\ \quad \text{\(\uparrow \)} \\ \quad \text{\(\uparrow \)} \\ \quad \text{\(\uparrow \)} \\ \quad \text{\(\uparrow \)} \\ \quad \text{\(\uparrow \)} \end{array}$$

machine (IM) in rotor field reference system operating at constant flux (Leonhard, 1996): Figure 2. Standard space state form characterization of the three phase Induction Motor. **Figure 2.** Standard space state form characterization of the three phase Induction Motor

#### **3. Optimal control problem statement**

The common objectives of the optimal control drive systems are: the smooth dynamic response, without oscillations; achieving the desired steady state; the fast compensation of the load torque; energy minimization. Taking into consideration the main objectives of the electrical drive system, the following Section describes the chosen quadratic performance criterion.

( ) <sup>1</sup> ( , , ,) () () () <sup>2</sup> *Hpxut tt t* = + +× ++ é ù ë û *T TT x (t) Qx(t) u (t)Ru(t) y (t) Ax Bu Gw* (4)

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

( ) \* *u t*( ) *t* - = - *1 T R By* (5)

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

345

)

*<sup>x</sup>* (7)

<sup>0</sup> <sup>0</sup> , and *q* is the angle of the rotor magnetizing flux.

<sup>o</sup> (6)

The condition **R>0** assures the inversability of the eigenvectors and that the optimal solution has a minimum (Gaiceanu, 2002). The optimal control minimizing the Hamiltonian (4) is given

By using the canonical system, the solutions of the costate *y***(***t***)** and state vectors *x***(***t***)** can be


*1 T T x(t) A BR B x(t) G w(t Q A y(t) 0 y(t)*

= = - é ù ë û *1 1 1 1 t t λ(t) y(t ) S x(t ) x*

<sup>0</sup> ,

The integration of the canonical system leads to the well-known matrix Riccati differential equation and the associate vectorial equation (Rosu, 1985), (Gaiceanu, 2002), (Athans, *et al.*, 2006). The integration of these two equations is a very difficult task because the Riccati equation is nonlinear, its solution is recursive (which can be calculated backward in time) (Gattami, *et al.*, 2005), (Jianqiang, *et al.*, 2007), (Rosu, 1985). Moreover, the backward computation needs to know a priori the variation of the perturbation vector **w**(*t*) during the control interval [0, *t*1].

in which *y* **(***t***)***∈ℜ<sup>2</sup>*

by

obtained:


*y*(*t*1)=*S x*(*t*1)− *x*<sup>1</sup> =100

is the costate vector.

o


*ω*(*t*1)

where the weighting matrix *<sup>S</sup>* <sup>=</sup> <sup>100</sup> <sup>0</sup>

The boundary conditions of the canonical system (5) are:

¶

*<sup>q</sup>*(*t*1) <sup>−</sup> 15707,96

In most electrical drives the last condition cannot be accomplished.

¶=

#### **3.1. The quadratic performance criterion**

The functional cost or the quadratic performance criterion is chosen such that the electrical input power, the energy expended in the stator windings of the electrical machines is mini‐ mized, and the cost of the control is weighted. Therefore, the performance functional quadratic criterion (Rosu, 1985), (Gaiceanu, 2002), (Athans, 2006), is as follows:

$$J = \frac{1}{2} \left[ \mathbf{x}(t\_1) - \mathbf{x}\_1 \right]^T \mathbf{S} \left[ \mathbf{x}(t\_1) - \mathbf{x}\_1 \right] + \frac{1}{2} \int\_{t\_0}^{t\_1} \left[ \mathbf{x}^T(t) \, \mathbf{Q} \mathbf{x}(t) + \mathbf{u}^T(t) \, \mathbf{R} \mathbf{u}(t) \right] dt \tag{3}$$

in which, taking into consideration a starting of the IM, the initial conditions consists of the initial time *t*0=0[s] and the initial state vector *x*(*t*0) <sup>=</sup> *<sup>ω</sup>*(*t*0) *<sup>q</sup>*(*t*0) *<sup>T</sup>* <sup>=</sup> <sup>0</sup> <sup>0</sup> *<sup>T</sup>* , *<sup>T</sup>*-being the transpose symbol of the initial state vector. The final condition is the required final state, **x**1. The components of the **x**<sup>1</sup> are: the desired angular velocity, *ω*\* m, and a free rotor magnetizing angle value, *q\** =0, i.e *x*<sup>1</sup> <sup>=</sup> *<sup>ω</sup><sup>m</sup>* \* 0 ; the final time, *t*1=0.9[s]. The value of the final time is established by the feasibility condition, i.e. the value of it is taken from the obtained final time of the starting process with conventional control (Proportional Integral controllers). In order to exist and to obtain a unique optimal solution, the weighting matrices must be chosen according to: **S**≥0, **R**>0, **Q** ≥0. The weighted matrix **S** minimizes the square error between the reached state *x*(*t*1), and the desired state, **x**1, in the fixed time *t*1, the weighted matrix **R** minimizes the control effort, *u*(*t*)=*i qs* (*t*), and the weighted matrix **Q** minimizes the expended energy in the motor windings. The optimal control problem taken into account in this chapter is with free-end point, fixed time and unconstrained. The restrictions of the magnitude for the control vector, and state vectors are managed through the design process, by an adequate choice of the weighting matrices.

#### **3.2. The solution of the optimal control problem.**

The existence and the unicity of the optimal control problem solution are based on the controllability and observability of the system (1). The system (1) being controllable and completely observable, the weighting matrices are chosen such that **Q, S ≥ 0** are positive semidefinite, and matrix **R>0** is positive definite (Gaiceanu, 2002). In this way, the existence and the unicity of the optimal control problem solution are assured.

By using variational method, the Hamiltonian of the optimal control problem is given by

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives http://dx.doi.org/10.5772/57521 345

$$H(p, \mathbf{x}, u, t) = \frac{1}{2} \left[ \mathbf{x}^T(\mathbf{t}) \, \mathbf{Q} \mathbf{x}(\mathbf{t}) + \mathbf{u}^T(\mathbf{t}) \mathbf{R} \mathbf{u}(\mathbf{t}) + \mathbf{y}^T(\mathbf{t}) \cdot \left( \mathbf{A} \mathbf{x}(t) + \mathbf{B} \mathbf{u}(t) + \mathbf{G} \mathbf{w}(t) \right) \right] \tag{4}$$

in which *y* **(***t***)***∈ℜ<sup>2</sup>* is the costate vector.

**3. Optimal control problem statement**

344 MATLAB Applications for the Practical Engineer

**3.1. The quadratic performance criterion**

=0, i.e *x*<sup>1</sup> <sup>=</sup> *<sup>ω</sup><sup>m</sup>*

angle value, *q\**

effort, *u*(*t*)=*i*

weighting matrices.

*qs*

The common objectives of the optimal control drive systems are: the smooth dynamic response, without oscillations; achieving the desired steady state; the fast compensation of the load torque; energy minimization. Taking into consideration the main objectives of the electrical drive system, the following Section describes the chosen quadratic performance criterion.

The functional cost or the quadratic performance criterion is chosen such that the electrical input power, the energy expended in the stator windings of the electrical machines is mini‐ mized, and the cost of the control is weighted. Therefore, the performance functional quadratic

1

in which, taking into consideration a starting of the IM, the initial conditions consists of the initial time *t*0=0[s] and the initial state vector *x*(*t*0) <sup>=</sup> *<sup>ω</sup>*(*t*0) *<sup>q</sup>*(*t*0) *<sup>T</sup>* <sup>=</sup> <sup>0</sup> <sup>0</sup> *<sup>T</sup>* , *<sup>T</sup>*-being the transpose symbol of the initial state vector. The final condition is the required final state, **x**1.

by the feasibility condition, i.e. the value of it is taken from the obtained final time of the starting process with conventional control (Proportional Integral controllers). In order to exist and to obtain a unique optimal solution, the weighting matrices must be chosen according to: **S**≥0, **R**>0, **Q** ≥0. The weighted matrix **S** minimizes the square error between the reached state *x*(*t*1), and the desired state, **x**1, in the fixed time *t*1, the weighted matrix **R** minimizes the control

windings. The optimal control problem taken into account in this chapter is with free-end point, fixed time and unconstrained. The restrictions of the magnitude for the control vector, and state vectors are managed through the design process, by an adequate choice of the

The existence and the unicity of the optimal control problem solution are based on the controllability and observability of the system (1). The system (1) being controllable and completely observable, the weighting matrices are chosen such that **Q, S ≥ 0** are positive semidefinite, and matrix **R>0** is positive definite (Gaiceanu, 2002). In this way, the existence

By using variational method, the Hamiltonian of the optimal control problem is given by

*<sup>t</sup> J t <sup>t</sup> t t t t dt* é ù = - -+ é ùé ù <sup>+</sup> ë ûë û ò ë û *x x S x x x Qx u Ru* (3)

(*t*), and the weighted matrix **Q** minimizes the expended energy in the motor

; the final time, *t*1=0.9[s]. The value of the final time is established

m, and a free rotor magnetizing

1 1 () () () () () () 2 2 *T t T T*

criterion (Rosu, 1985), (Gaiceanu, 2002), (Athans, 2006), is as follows:

<sup>0</sup> 11 11

The components of the **x**<sup>1</sup> are: the desired angular velocity, *ω*\*

\*

0

**3.2. The solution of the optimal control problem.**

and the unicity of the optimal control problem solution are assured.

The condition **R>0** assures the inversability of the eigenvectors and that the optimal solution has a minimum (Gaiceanu, 2002). The optimal control minimizing the Hamiltonian (4) is given by

$$
\mu^\*(t) = -\mathbb{R}^{-1}\mathbb{B}^T\mathcal{Y}(t)\tag{5}
$$

By using the canonical system, the solutions of the costate *y***(***t***)** and state vectors *x***(***t***)** can be obtained:

$$
\begin{bmatrix} \stackrel{\circ}{\mathbf{x}}(\mathbf{t}) \\ \stackrel{\circ}{\mathbf{y}}(\mathbf{t}) \end{bmatrix} = \begin{bmatrix} A & -\mathbf{B}\mathbf{R}^{-1}\mathbf{B}^{T} \\ -\mathbf{Q} & -A^{T} \end{bmatrix} \begin{bmatrix} \mathbf{x}(\mathbf{t}) \\ \mathbf{y}(\mathbf{t}) \end{bmatrix} + \begin{bmatrix} \mathbf{G} \\ \mathbf{0} \end{bmatrix} w(\mathbf{t}) \tag{6}
$$

The boundary conditions of the canonical system (5) are:



$$\mathbf{y}(\mathbf{t}\_1) = \frac{\partial \lambda(\mathbf{t})}{\partial \mathbf{x}}\bigg|\_{\mathbf{t}=\mathbf{t}\_1} = \mathbf{S}\Big[\mathbf{x}(\mathbf{t}\_1) - \mathbf{x}\_1\Big] \tag{7}$$

$$\mathbf{y}(t\_1) = \mathbf{S} \begin{bmatrix} \mathbf{x}(t\_1) - \mathbf{x}\_1 \end{bmatrix} = 100 \begin{bmatrix} \omega(t\_1) \\ q(t\_1) \end{bmatrix} \quad - \quad \begin{bmatrix} 15707,96 \\ 0 \end{bmatrix} \mathbf{y}$$

where the weighting matrix *<sup>S</sup>* <sup>=</sup> <sup>100</sup> <sup>0</sup> <sup>0</sup> <sup>0</sup> , and *q* is the angle of the rotor magnetizing flux.

The integration of the canonical system leads to the well-known matrix Riccati differential equation and the associate vectorial equation (Rosu, 1985), (Gaiceanu, 2002), (Athans, *et al.*, 2006). The integration of these two equations is a very difficult task because the Riccati equation is nonlinear, its solution is recursive (which can be calculated backward in time) (Gattami, *et al.*, 2005), (Jianqiang, *et al.*, 2007), (Rosu, 1985). Moreover, the backward computation needs to know a priori the variation of the perturbation vector **w**(*t*) during the control interval [0, *t*1]. In most electrical drives the last condition cannot be accomplished.

A nonrecursive solution of the Riccati equation has been developed by using two linear transformations (Rosu, *et al.*, 1998). *The first transformation* changes actual time *t* into *t1-t,* that is time remaining until the end of the optimal process.

$$
\tau = t\_1 - t = 0.9 - t\_\prime \tag{8}
$$

In the Appendix 3 the proper position of the eigenvalues and the eigenvectors of the linear system (1) is provided. By using the program sequence from Appendix 3, the stability of the

> 0.9877 -0.5906 -0.9740 -0.5702 -0.1567 0.8070 -0.2268 -0.8215 0.0004 -0.0001 0.0004 0.0000 -0.0000 0.0011 0.0000 0.0011

*W* (14)

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

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

347

 0.0005 0.0001 1.3542 0.3153 0.0000 0.0006 0.3158 0.4549

é ù ê ú

0.0005 0.0000 1.3732 0.2178 0.0000 0.0006 0.3271 0.4469

*W* (15)

*n q* (16)

*q n* (17)

t

<sup>o</sup> (18)

t

<sup>o</sup> (19)


é ù ê ú

ë û


<sup>1</sup> () () , () ()

() () () ()

éù é ù êú ê ú = ëû ë û *p m <sup>W</sup>*


 t

 t

 t

 t

By introducing (15) and (16) transformations in (8), the following differential system is

() () ( ) ( ) ( )

*n 0*

t

1 1 ( ) ( ) ( ) ( ) ( )


 t

t

*W MW W r n 0*

ê ú é ù éù <sup>=</sup> ê ú êú - ê ú ë û ëû ê ú ë û *m mG W MW r*

Therefore, the canonical system can be written in terms of **n**(*τ*), and **m**(*τ*) vectors:

ê ú é ù éù <sup>=</sup> ê ú êú - ê ú ë û ëû ê ú ë û *m mG*

 t

system is assured and the corresponding matrix of the eigenvectors is obtained:

=

By knowing the inverse of the eigenvector matrix **W**-1

1 1.0 003 \*

the new vectors, **m** (*τ*), **n**(*τ*), **p**(*τ*), **q**(*τ*) are introduced as follows:

t

é ù

o

t

*n*

t

é ù

o

t

*n*

t

t

t

t

*e*

obtained:

In order to use only the negative eigenvalues for the computation of the fundamental matrix of the system (1), *the second transformation* is made. Therefore, by defining new vectors **p**(*τ*), **q**(*τ*), **r**(*τ*),

$$
\begin{bmatrix}
\stackrel{\circ}{\mathcal{p}}(\tau) \\
\stackrel{\circ}{\mathcal{q}}(\tau)
\end{bmatrix} = \begin{bmatrix}
\mathbf{Q} & A^{T}
\end{bmatrix} \begin{bmatrix}
\mathcal{p}(\tau) \\
\mathcal{q}(\tau)
\end{bmatrix} - \begin{bmatrix}
\mathbf{G} \\
\mathbf{0}
\end{bmatrix} r(\tau) \tag{9}
$$

and by noting the canonical matrix (4×4) of the system (1) as:

$$
\begin{bmatrix}
\mathbf{Q} & \mathbf{A}^{T}
\end{bmatrix} = \mathbf{M} \tag{10}
$$

the eigenvalues *λ<sup>i</sup>* , *i* =1 ¯ , 4 of the canonical matrix **M** can be computed with:

$$\det\begin{bmatrix} \lambda I - \mathbf{M} \end{bmatrix} = 0\tag{11}$$

By using the developed Matlab file from Appendix 3, the appropriate negative eigenvalues position placed on the main diagonal of the eigenvalues matrix, **D**, with diag(*D*)={*λ*1, *λ*2, *λ*3, *λ*4}, can be obtained:

$$
\begin{bmatrix}
\lambda\_1 & 0 & 0 & 0 \\
0 & \lambda\_2 & 0 & 0 \\
0 & 0 & -\lambda\_1 & 0 \\
0 & 0 & 0 & -\lambda\_2
\end{bmatrix} = \begin{bmatrix}
5.1663 & 0 & 0 & 0 \\
0 & 0.7205 & 0 & 0 \\
0 & 0 & -5.1663 & 0 \\
0 & 0 & 0 & -0.7205
\end{bmatrix} = \begin{bmatrix}
\Lambda & \mathbf{0} \\
\mathbf{0} & -\Lambda
\end{bmatrix} \tag{12}
$$

in which

$$
\Lambda = \begin{bmatrix} \lambda\_1 & 0 \\ 0 & \lambda\_2 \end{bmatrix} = \begin{bmatrix} 5.1663 & 0 \\ 0 & 0.7205 \end{bmatrix} . \tag{13}
$$

In the Appendix 3 the proper position of the eigenvalues and the eigenvectors of the linear system (1) is provided. By using the program sequence from Appendix 3, the stability of the system is assured and the corresponding matrix of the eigenvectors is obtained:

$$\mathbf{W} = \begin{bmatrix} 0.9877 & -0.5906 & -0.9740 & -0.5702 \\ -0.1567 & 0.8070 & -0.2268 & -0.8215 \\ 0.0004 & -0.0001 & 0.0004 & 0.0000 \\ -0.0000 & 0.0011 & 0.0000 & 0.0011 \end{bmatrix} \tag{14}$$

By knowing the inverse of the eigenvector matrix **W**-1

A nonrecursive solution of the Riccati equation has been developed by using two linear transformations (Rosu, *et al.*, 1998). *The first transformation* changes actual time *t* into *t1-t,* that

In order to use only the negative eigenvalues for the computation of the fundamental matrix of the system (1), *the second transformation* is made. Therefore, by defining new vectors **p**(*τ*),

> <sup>1</sup> ( ) ( ) ( ) ( ) ( ) *T T*

> > 1 *T T* - é ù ê ú <sup>=</sup> ê ú ë û *A BR B*

deté ù l

00 0 5.1663 0 0 0

00 0 0 0 0 -0.7205

2

l


1

l

0 00 0 0.7205 0 0 , 00 0 0 0 -5.1663 0

<sup>0</sup> 5.1663 0 . <sup>0</sup> 0 0.7205

é ù é ù L= = ê ú ê ú ë û ë û (13)

é ù é ù ê ú ê ú é ù <sup>L</sup> = = ê ú - ë û -L

By using the developed Matlab file from Appendix 3, the appropriate negative eigenvalues position placed on the main diagonal of the eigenvalues matrix, **D**, with


 t

t

*M*

, 4 of the canonical matrix **M** can be computed with:

= -= - *tt t* 0.9 , (8)

t

*Q A* (10)


*0*

*<sup>0</sup>* (12)

<sup>o</sup> (9)

1 t

is time remaining until the end of the optimal process.

346 MATLAB Applications for the Practical Engineer

t

t

and by noting the canonical matrix (4×4) of the system (1) as:

*q*

, *i* =1

diag(*D*)={*λ*1, *λ*2, *λ*3, *λ*4}, can be obtained:

1

l

2

l

¯

o

**q**(*τ*), **r**(*τ*),

the eigenvalues *λ<sup>i</sup>*

1

l

in which

2

l

$$\begin{array}{cccc} \textbf{W-1} = \textbf{1.0e} + \textbf{0.03} & \textbf{4.0001} & \textbf{1.3542} & \textbf{0.3153} \\ \textbf{0.0000} & \textbf{0.0006} & \textbf{0.3158} & \textbf{0.4549} \\ \textbf{-0.0005} & \textbf{0.0000} & \textbf{1.3732} & \textbf{-0.2178} \\ \textbf{0.0000} & \textbf{-0.0006} & \textbf{-0.3271} & \textbf{0.4469} \end{array} \tag{15}$$

the new vectors, **m** (*τ*), **n**(*τ*), **p**(*τ*), **q**(*τ*) are introduced as follows:

$$
\begin{bmatrix} \mathfrak{m}(\tau) \\ \mathfrak{m}(\tau) \end{bmatrix} = \mathsf{W}^{-1} \begin{bmatrix} \mathfrak{p}(\tau) \\ \mathfrak{q}(\tau) \end{bmatrix}' \tag{16}
$$

$$
\begin{bmatrix} p(\boldsymbol{\tau}) \\ q(\boldsymbol{\tau}) \end{bmatrix} = \mathbf{W} \begin{bmatrix} m(\boldsymbol{\tau}) \\ n(\boldsymbol{\tau}) \end{bmatrix} \tag{17}
$$

By introducing (15) and (16) transformations in (8), the following differential system is obtained:

$$\mathbf{W} \begin{bmatrix} \stackrel{\circ}{m}(\tau) \\ \stackrel{\circ}{m}(\tau) \\ \stackrel{\circ}{n(\tau)} \end{bmatrix} = \mathbf{M} \mathbf{W} \begin{bmatrix} m(\tau) \\ n(\tau) \end{bmatrix} \quad - \quad \begin{bmatrix} G \\ \mathbf{0} \end{bmatrix} r(\tau) \tag{18}$$

Therefore, the canonical system can be written in terms of **n**(*τ*), and **m**(*τ*) vectors:

$$
\begin{bmatrix} \stackrel{\circ}{m}(\tau) \\ \stackrel{\circ}{m}(\tau) \\ \stackrel{\circ}{n(\tau)} \end{bmatrix} = \mathbf{W}^{-1} \mathbf{M} \mathbf{W} \begin{bmatrix} m(\tau) \\ n(\tau) \end{bmatrix} - \mathbf{W}^{-1} \begin{bmatrix} G \\ \mathbf{0} \end{bmatrix} r(\tau) \tag{19}
$$

By using the well-known relation

$$\mathbf{W}^{-1}\mathbf{M}\mathbf{W} = \begin{bmatrix} \Lambda & 0\\ 0 & -\Lambda \end{bmatrix} \tag{20}$$

where **I** is the identity matrix.

t

t

t

t

é - ê

( )

5.1663

t

where

where,

t

×

The obtained solution of the canonical system (23) is:

t

×

( )

1 e 00 0

t

×

0 1e 0 0 0 0 1e 0 0 0 0 1e

<sup>1</sup> 1 0.7205 2 2 5.1663 <sup>1</sup> <sup>1</sup> 2 0.7205 2

*m m m m n n n n*

t

×

e 000 ( ) (0) ( ) 0e 0 0 (0) ( ) (0) 0 0e 0 ( ) (0) 0 0 0e

é ù é ù é ù ê ú ê ú ê ú ê ú ê ú <sup>=</sup> × + ê ú ê ú ê ú ê ú ê ú ë û ê ú ë û ë û

( )



the first equation of the system **(**24**)** becomes:

t 0.7205

( )


t

( )

<sup>ù</sup> é ù <sup>ú</sup> ê ú ê ú ê ú × × ê ú



(0) 0 () 0


t

2

The transversality condition (6), in terms of **m** and **n** vectors, is as follows:

1 12 22 21 11

é ù é ù é ù é ù - é ù = ×+ ê ú ê ú× × ê ú ê ú ê ú ë û ê ú ë û ë û - ë û ë û *me m I e H*

( ) 0 (0) 0

t 5.1663

t

The negative exponentials can be obtained by extracting **m**(0) function on the **m**(τ). In this way,


1 0 0 1 é ù <sup>=</sup> ê ú ë û

The negative exponentials of the system **(**25**)** guarantee the stability of the system (1).


( )

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives


t

( )

2 1

t

*n n <sup>e</sup> I e <sup>H</sup>* (26)

2 2

*I* (27)

(0) (0) *<sup>f</sup> n Em Fx* =× +× (28)

 1.0215 -0.0320 -0.0127 -0.9811

*E SW W W SW* (29)

t

0.7205


t


23.0137

ê ú ê ú ë û <sup>û</sup>


*r*

*r*

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

(25)

349

( )

5.1663

and by noting the new vector **H** as:

$$\mathbf{W}^{-1}\begin{bmatrix}\mathbf{G}\\\mathbf{0}\end{bmatrix} = \mathbf{H} = \begin{bmatrix} -52.6495\\ -0.2471\\ -47.8790\\ 23.0137 \end{bmatrix} \tag{21}$$

the new form of the canonical system can be obtained:

$$
\begin{bmatrix} \stackrel{\circ}{m}(\tau) \\ \stackrel{\circ}{m}(\tau) \end{bmatrix} = \begin{bmatrix} \Lambda & 0 \\ 0 & -\Lambda \end{bmatrix} \begin{bmatrix} m(\tau) \\ m(\tau) \end{bmatrix} - H \cdot r(\tau) \tag{22}
$$

The solution of the new canonical system (21) can be obtained by using the following equation:

$$
\begin{bmatrix} \mathfrak{m}(\tau) \\ \mathfrak{m}(\tau) \end{bmatrix} = \begin{bmatrix} e^{\Lambda \tau} & 0 \\ 0 & e^{-\Lambda \tau} \end{bmatrix} \begin{bmatrix} \mathfrak{m}(0) \\ \mathfrak{m}(0) \end{bmatrix} - \begin{bmatrix} e^{\Lambda \beta} & 0 \\ 0 & e^{-\Lambda \beta} \end{bmatrix} H \mathbf{r}(\beta) \mathbf{d}\beta \tag{23}
$$

As it can be observed the solution of the canonical system (22) solution contains both the positive, *e Λτ* , and the negative exponentials, *e* <sup>−</sup>*Λτ* . It is well-known that the positive exponen‐ tials produce the instability of the solution. Moreover, in order to obtain the analytical solution of the system (22) the load torque (*r*(τ) or **w**(*t*)) at the final time, *t*1, must be known in advance. For the most used electrical drive system this condition is hardly accomplished because the load torque cannot be known in advance, at the final time.

A numerical solution of the optimal problem is proposed. The advantage of maintaining the sampled signal at the same value during the sampled period is that the zero order hold (ZOH) offers a brilliant solution to know the load torque at the final time. Therefore, under this assumption, the solution **(**22**)** is computed each sampled time, *T,* and the perturbation vector **r** has constant value during it:

$$
\begin{bmatrix} m(\tau) \\ m(\tau) \end{bmatrix} = \begin{bmatrix} e^{\Lambda \tau} & 0 \\ 0 & e^{-\Lambda \tau} \end{bmatrix} \begin{bmatrix} m(0) \\ m(0) \end{bmatrix} + \begin{bmatrix} I - e^{\Lambda \tau} & 0 \\ 0 & I - e^{-\Lambda \tau} \end{bmatrix} \begin{bmatrix} H\_1 \\ H\_2 \end{bmatrix} r \tag{24}
$$

where **I** is the identity matrix.

By using the well-known relation

348 MATLAB Applications for the Practical Engineer

and by noting the new vector **H** as:

<sup>1</sup> 0 0 - é ù <sup>L</sup> <sup>=</sup> ê ú

ê ú é ù ê ú ê ú = = ê ú ë û ê ú

() 0 () ( ) 0 () ( )

ê ú é ùé ù <sup>L</sup> <sup>=</sup> - × ê ú ê úê ú ë ûë û -L ê ú

*m m*

t

*n n e e*

L L

éù é ù éù éù êú êú = êú ê ú ëû ëû ëû ë û <sup>ò</sup> *me m e*

( ) 0 (0) 0

L L

( ) 0 0 (0)

t

 t

t

*n*

The solution of the new canonical system (21) can be obtained by using the following equation:

0 ( ) 0 (0) <sup>0</sup> ( )d ( ) 0 0 (0)

As it can be observed the solution of the canonical system (22) solution contains both the

tials produce the instability of the solution. Moreover, in order to obtain the analytical solution of the system (22) the load torque (*r*(τ) or **w**(*t*)) at the final time, *t*1, must be known in advance. For the most used electrical drive system this condition is hardly accomplished because the

A numerical solution of the optimal problem is proposed. The advantage of maintaining the sampled signal at the same value during the sampled period is that the zero order hold (ZOH) offers a brilliant solution to know the load torque at the final time. Therefore, under this assumption, the solution **(**22**)** is computed each sampled time, *T,* and the perturbation vector

> -L -L éù éù é ùé ù - é ù êú êú = + ê úê ú ê ú ëû ëû ë ûë û - ë û *m e m Ie H*

 t

*n n <sup>e</sup> I e <sup>H</sup>* (24)


t b t

 b b b

. It is well-known that the positive exponen‐

1 2

 t *r*

(23)

*Hr*

<sup>o</sup> (22)

*H r*

1

*G W H*


the new form of the canonical system can be obtained:

t

é ù

o

t

ë û

t

, and the negative exponentials, *e* <sup>−</sup>*Λτ*

load torque cannot be known in advance, at the final time.

t

t

t

positive, *e Λτ*

**r** has constant value during it:

t

t

*n*

ë û -L


é ù

23.0137

ê ú ë û


*W MW* (20)

*<sup>0</sup>* (21)

The obtained solution of the canonical system (23) is:

$$
\begin{bmatrix} m\_1(\mathbf{r}) \\ m\_2(\mathbf{r}) \\ n\_1(\mathbf{r}) \\ n\_2(\mathbf{r}) \\ \end{bmatrix} = \begin{bmatrix} \mathbf{e}^{(5.1663)\tau} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{e}^{(0.7205)\tau} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{e}^{-(5.1663)\tau} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{e}^{-(0.7205)\tau} \end{bmatrix} \cdot \begin{bmatrix} m\_1(\mathbf{0}) \\ m\_2(\mathbf{0}) \\ n\_1(\mathbf{0}) \\ n\_2(\mathbf{0}) \end{bmatrix} + \\ \begin{bmatrix} 1 - \mathbf{e}^{(5.1663)\tau} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & 1 - \mathbf{e}^{(0.7205)\tau} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & 1 - \mathbf{e}^{(-5.1663)\tau} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{1} - \mathbf{e}^{(-5.1663)\tau} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & 1 - \mathbf{e}^{(0.7205)\tau} \end{bmatrix} \cdot \begin{bmatrix} -52.6495 \\ -0.2471 \\ -47.8790 \\ 23.0137 \end{bmatrix} \cdot r. \end{bmatrix} \tag{25}
$$

The negative exponentials can be obtained by extracting **m**(0) function on the **m**(τ). In this way, the first equation of the system **(**24**)** becomes:

$$
\begin{bmatrix} m(0) \\ m(\tau) \end{bmatrix} = \begin{bmatrix} e^{-\Lambda \cdot \tau} & 0 \\ 0 & e^{-\Lambda \cdot \tau} \end{bmatrix} \cdot \begin{bmatrix} m(\tau) \\ n(0) \end{bmatrix} + \begin{bmatrix} I\_2 - e^{-\Lambda \cdot \tau} & 0 \\ 0 & I\_2 - e^{-\Lambda \cdot \tau} \end{bmatrix} \cdot \begin{bmatrix} H\_1 \\ H\_2 \end{bmatrix} \cdot r \tag{26}
$$

where

$$I\_2 = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \tag{27}$$

The negative exponentials of the system **(**25**)** guarantee the stability of the system (1).

The transversality condition (6), in terms of **m** and **n** vectors, is as follows:

$$
\mathfrak{m}(0) = \mathbf{E} \cdot \mathfrak{m}(0) + \mathbf{F} \cdot \mathbf{x}\_f \tag{28}
$$

where,

$$E = \left[\mathbf{S}\mathbf{W}\_{12} - \mathbf{W}\_{22}\right]^{-1} \left[\mathbf{W}\_{21} - \mathbf{S}\mathbf{W}\_{11}\right] = \begin{bmatrix} 1.0215 & -0.0320 \\ -0.0127 & -0.9811 \end{bmatrix} \tag{29}$$

$$F = \left[\mathbf{S}\mathbf{W}\_{12} - \mathbf{W}\_{22}\right]^{-1}\mathbf{S} = \begin{bmatrix} -1.0504 & 0\\ 0.0404 & 0 \end{bmatrix} \tag{30}$$

By expressing **n**(*τ*) function on **m**(*τ*), the following equation can be deducted :

$$\mathfrak{m}(\tau) = \mathbf{Z}(\tau) \cdot \mathfrak{m}(\tau) + e^{-\Lambda \cdot \tau} \cdot \mathbf{E} \cdot \left( I\_2 - e^{-\Lambda \cdot \tau} \right) \cdot H\_1 \cdot r + e^{-\Lambda \cdot \tau} \cdot \mathbf{F} \cdot \mathbf{x}\_1 + \left( I\_2 - e^{-\Lambda \cdot \tau} \right) \cdot H\_2 \cdot r \tag{31}$$

where

$$\mathbf{Z}(\tau) = \mathbf{e}^{-\Lambda \cdot \tau} \cdot \mathbf{E} \cdot \mathbf{e}^{-\Lambda \cdot \tau} \tag{32}$$

The optimal control, at any moment *t*, is given by (Gaiceanu, 2002), (Athans, *et al.*, 2006)

Stator voltage eqs. d/q u\* q1 u\* A1 u\* B1 u\* C1

Power Inverter SVPWM

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

3~

351

iA1 iB1 iC1

> M 3~

IM mathematical model

u\* d1

x(t) q(t)

Load torque estimator 2/3 e -jq

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

ids(t) iqs(t)

TL(t) <sup>w</sup>(t)

sq) implementation into three phase IM Drive System supplied by the Voltage Power

w(t)

in which **P(t1-t)** is the solution of the matrix Riccati differential equation, MRDE, and the

The numerical solution supposes the knowledge of the disturbance **w**(*t*)=*TL*(*t*), which could be available by using torque estimator, as in (Rosu, *et al.*, 1998), or with torque sensor (more

In order to implement the optimal control system two electrical drive systems are combined: the conventional vector control and the optimal control drive. Both electrical drive systems operate in parallel. The activation of one of them depends on the speed error: *if* the speed error is zero (steady state operating regime), *then* the *conventional vector control* is active, *else* the

Optimal control implementation consists of loading nine software modules developed by the author: nameplate motor data (Fig.1), determination of the motor parameters (Appendix 1), determination of the six sectors vector components to be used in Space Vector Modulation

<sup>1</sup> 11 1 2 1 () ( ) () ( ) ( ) () *T TT ut t t t t t tt t* - -- =- - + - + - *R BP x R BK x R BK w* (37)

\* 1 1 1

expensive solution that can produce the faults of the electrical drive system).

matrices K1 and K2 are calculated via P(t1-t) (Rosu, *et al.*, 1998).


i \* qs(t)

i \* <sup>d</sup>s(t)

K2 (t1 -t)

\*

**4. Simulation and experimental results**

x1 R-1 BT

K1 (t1 -t)

**Figure 3.** Optimal control (*i*

Source Inverter

*optimal control* drive is active.

**4.1. Optimal control implementation**

Going back through the reverse transformations from **m**(*τ*), **n**(*τ*) vectors to **p**(*τ*), **q**(*τ*) vectors, from the used *τ* time to the current time *t*, the state **x**(*t***)** and costate **y**(*t*) expressions can be obtained. The costate vector at current time is as follows:

$$\mathbf{y}(t) = \mathbf{P}(t\_1 - t) \cdot \mathbf{x}(t) - \mathbf{K}\_1(t\_1 - t) \cdot \mathbf{x}\_1 - \mathbf{K}\_2(t\_1 - t) \cdot \mathbf{w} \tag{33}$$

where the matrix **P**(*t*1−*t*) is the solution of the matrix Riccati diferential equation, and the matrices**K**1(*t*1−*t*) and **K**2(*t*1*-t***)** vector have the form:

$$\mathbf{K}\_1(t\_1 - t) = -\boldsymbol{\sigma}(t\_1 - t) \cdot \mathbf{e}^{-\boldsymbol{\Lambda}\cdot\left\{t\_1 - t\right\}} \cdot \mathbf{F}\_\prime \tag{34}$$

$$\mathbf{K}\_2(t\_1 - t) = -\boldsymbol{\sigma}(t\_1 - t) \cdot \left[ \mathbf{e}^{-\boldsymbol{\Lambda}\cdot\left(t\_1 - t\right)} \cdot \mathbf{E} \cdot \left( \mathbf{I} - \mathbf{e}^{-\boldsymbol{\Lambda}\cdot\left(t\_1 - t\right)} \right) \cdot \mathbf{H}\_1 + \left( \mathbf{I} - \mathbf{e}^{-\boldsymbol{\Lambda}\cdot\left(t\_1 - t\right)} \right) \cdot \mathbf{H}\_2 \right],\tag{35}$$

The **v** matrix is as follows:

$$\mathbf{w}(\tau) = \left[\mathbf{W}\_{22} - \mathbf{P}(\tau)\mathbf{W}\_{12}\right]. \tag{36}$$

The solution of the optimal control problem has three components (Fig.3):


Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives http://dx.doi.org/10.5772/57521 351

**Figure 3.** Optimal control (*i* \* sq) implementation into three phase IM Drive System supplied by the Voltage Power Source Inverter

The optimal control, at any moment *t*, is given by (Gaiceanu, 2002), (Athans, *et al.*, 2006)

$$\boldsymbol{\mu}^\*(t) = -\mathbf{R}^{-1}\mathbf{B}^T\mathbf{P}(t\_1 - t)\mathbf{x}(t) + \mathbf{R}^{-1}\mathbf{B}^T\mathbf{K}\_1(t\_1 - t)\mathbf{x}\_1 + \mathbf{R}^{-1}\mathbf{B}^T\mathbf{K}\_2(t\_1 - t)\mathbf{w}(t) \tag{37}$$

in which **P(t1-t)** is the solution of the matrix Riccati differential equation, MRDE, and the matrices K1 and K2 are calculated via P(t1-t) (Rosu, *et al.*, 1998).

The numerical solution supposes the knowledge of the disturbance **w**(*t*)=*TL*(*t*), which could be available by using torque estimator, as in (Rosu, *et al.*, 1998), or with torque sensor (more expensive solution that can produce the faults of the electrical drive system).

#### **4. Simulation and experimental results**

In order to implement the optimal control system two electrical drive systems are combined: the conventional vector control and the optimal control drive. Both electrical drive systems operate in parallel. The activation of one of them depends on the speed error: *if* the speed error is zero (steady state operating regime), *then* the *conventional vector control* is active, *else* the *optimal control* drive is active.

#### **4.1. Optimal control implementation**

1


() () () ( 2 1 ) *r r* 12 2 ( )

t

Going back through the reverse transformations from **m**(*τ*), **n**(*τ*) vectors to **p**(*τ*), **q**(*τ*) vectors, from the used *τ* time to the current time *t*, the state **x**(*t***)** and costate **y**(*t*) expressions can be

where the matrix **P**(*t*1−*t*) is the solution of the matrix Riccati diferential equation, and the

<sup>22</sup> <sup>12</sup> () () .

**1.** the state feedback *<sup>R</sup>* <sup>−</sup><sup>1</sup> <sup>⋅</sup> *<sup>B</sup><sup>T</sup> <sup>P</sup>*(*t*<sup>1</sup> <sup>−</sup>*t*) <sup>⋅</sup> *<sup>x</sup>*(*t*), in order to assure the stability of the system;

**2.** the compensating feedforward load torque *<sup>R</sup>* <sup>−</sup><sup>1</sup> <sup>⋅</sup> *<sup>B</sup><sup>T</sup> <sup>K</sup>*2(*t*<sup>1</sup> <sup>−</sup>*t*) <sup>⋅</sup>*w*, for fast compensation of

**3.** the reference component *<sup>R</sup>* <sup>−</sup><sup>1</sup> <sup>⋅</sup> *<sup>B</sup><sup>T</sup> <sup>K</sup>*1(*t*<sup>1</sup> <sup>−</sup>*t*) <sup>⋅</sup> *<sup>x</sup>*1, to track accurately the reference signal.

 t= - é ù

t

The solution of the optimal control problem has three components (Fig.3):


*F SW W S* (30)


 t

ë û

 t


 t

<sup>1</sup> 11 1 21 *y P x K xK w* () ( ) () ( ) ( ) *t tt t tt tt* = -× - -× - -× (33)

1 1 <sup>1</sup> ()() , *t t tt tt* -L× - *K ve F* - =- - × × (34)

ë û *v W PW* (36)

( ) <sup>1</sup>

( ) ( ) ( ) ( ) ( ) 11 1 2 1 <sup>1</sup> <sup>1</sup> <sup>2</sup> ()() , *t t t t t t tt tt* é ù -L× - -L× - -L× - - =- - × × × - × + - × ê ú ë û *K v e EIe H Ie H* (35)

12 22

tt

( )

t

obtained. The costate vector at current time is as follows:

matrices**K**1(*t*1−*t*) and **K**2(*t*1*-t***)** vector have the form:

the main perturbation, i.e. the load torque;

t

where

 tt

350 MATLAB Applications for the Practical Engineer

The **v** matrix is as follows:

By expressing **n**(*τ*) function on **m**(*τ*), the following equation can be deducted :

Optimal control implementation consists of loading nine software modules developed by the author: nameplate motor data (Fig.1), determination of the motor parameters (Appendix 1), determination of the six sectors vector components to be used in Space Vector Modulation (SVM) (Appendix 2), standard space state form characterization of the three phase Induction Motor (Fig.2), the proper arrangement of the eigenvalues and eigenvectors (Appendix 3), optimal control implementation in MATLAB (Appendix 4), the main program (Appendix 5), graphical representation (Appendix 6), Space Vector Modulation –Matlab function **svm\_mod** (Appendix 7). By using the Matlab code from Appendix 8 the graphical representation of the simulation results can be obtained (Figures 4-7).

0 0.2 0.4 0.6 0.8

0 0.2 0.4 0.6 0.8

**Figure 5.** Optimal control, speed, longitudinal and transversal voltage components

time[s]

time[s]

usd=f(t)

0 0.2 0.4 0.6 0.8

0 0.2 0.4 0.6 0.8

time[s]

time[s]

usq=f(t)

speed=f(t)

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

353

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

usq [V]

voltage loci loculus


speed [rot/min]

optimal control isq=f(t)

0


usd [V]


**Figure 6.** Voltage loci loculus for a starting process




0

100

200

300

400

2

4

isq [A]

6

**Figure 4.** Three phase SVM supply voltages

#### **Remark**

The SVM increases the DC-link voltage usage, the maximum output voltage based on the space vector theory is 2 / 3 times higher than conventional sinusoidal modulation (Pulse Width Modulation). The purpose of the SVPWM technique is to approximate the reference voltage vector *Uout* by a combination of the eight switching patterns. The implemented Matlab sequence for obtaining the pulse pattern of the optimal SVM modulator (Appendix 2, Appen‐ dix 7) shown above is based on the developed methodology by the author (Gaiceanu, 2002). In optimal control case, there is no speed overshoot, and smoothly state response (speed and angle of the rotor magnetizing current) have been obtained (Fig.5, Fig.7).

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives http://dx.doi.org/10.5772/57521 353

**Figure 5.** Optimal control, speed, longitudinal and transversal voltage components

(SVM) (Appendix 2), standard space state form characterization of the three phase Induction Motor (Fig.2), the proper arrangement of the eigenvalues and eigenvectors (Appendix 3), optimal control implementation in MATLAB (Appendix 4), the main program (Appendix 5), graphical representation (Appendix 6), Space Vector Modulation –Matlab function **svm\_mod** (Appendix 7). By using the Matlab code from Appendix 8 the graphical representation of the

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

stator voltage of phase A [V]

time[s]

stator voltage of phase B [V]

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

time[s]

stator voltage of phase C [V]

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

time[s]

The SVM increases the DC-link voltage usage, the maximum output voltage based on the space vector theory is 2 / 3 times higher than conventional sinusoidal modulation (Pulse Width Modulation). The purpose of the SVPWM technique is to approximate the reference voltage vector *Uout* by a combination of the eight switching patterns. The implemented Matlab sequence for obtaining the pulse pattern of the optimal SVM modulator (Appendix 2, Appen‐ dix 7) shown above is based on the developed methodology by the author (Gaiceanu, 2002). In optimal control case, there is no speed overshoot, and smoothly state response (speed and

angle of the rotor magnetizing current) have been obtained (Fig.5, Fig.7).

simulation results can be obtained (Figures 4-7).


352 MATLAB Applications for the Practical Engineer



**Figure 4.** Three phase SVM supply voltages

**Remark**

usA=f(t)

usB=f(t)

usC=f(t)

**Figure 6.** Voltage loci loculus for a starting process

**Figure 7.** Optimal control *iqs*=*f*(*t*), current *imr*(*t*), load torque *TL*(*t*)(elevator), rotor magnetizing angle of the rotor field *q*(*t*) (left to right)

**Figure 8.** Simulink implementation of the rotor field orientated vector control of the three-phase induction motor

[-2\*In, 2\*In]

P+I/s

P+I/s

PI controller

PI controller1

**Figure 9.** Simulink implementation of the speed and flux Proportional-Integral (PI) regulators

Saturation

Saturation1

2 ID\* iq\*

COS(theta)

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

355

SIN(theta)

FOC mechanism

f lux

om

theta2

theta1

1 IQ\*

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

omref\*

3 flux\_ref

1 om\_ref

> 1 Ts.s+1 Sample&Hold1

4 flux

2 om

1 Ts.s+1 Sample&Hold

#### **4.2. Conventional vector control**

The conventional vector control consists of the rotor flux orientated drive system. The vector control is implemented in Simulink as following. The control system is in cascaded manner and consists of two major loops: the speed loop and the magnetic flux loop. The torque control loop is the minor loop of the speed control one. The references of the control loop are: the desired angular velocity for the speed loop, *ω*\* , the imposed stator current torque component for the torque loop, *i* \* *qs*, and the desired stator current flux component for the magnetization flux loop, *i* \* *ds*.

#### **Tuning of the speed and flux controllers**

By knowing the nameplate motor data (Fig.1), the tuning of the PI speed and flux controllers parameters can be done (see the developed Matlab file shown in the Appendix 8). The conventional vector control supposes the use of the two independent control axes: *d*-axis, along with the magnetizing flux direction, and *q* axis, in quadrature with *d*-axis. The most used field orientated control is that of the rotor orientated magnetizing flux (Fig.8). The magnetization flux and the torque can be controlled independently.

0 0.2 0.4 0.6 0.8

0 0.2 0.4 0.6 0.8

time[s]

time[s]

TL=f(t)

0 0.2 0.4 0.6 0.8

0 0.2 0.4 0.6 0.8

time[s]

, the imposed stator current torque component

*qs*, and the desired stator current flux component for the magnetization

time[s]

angle of the rotor magnetizing flux q=f(t)

magnetizing current=f(t)

0

q

**Figure 7.** Optimal control *iqs*=*f*(*t*), current *imr*(*t*), load torque *TL*(*t*)(elevator), rotor magnetizing angle of the rotor field

The conventional vector control consists of the rotor flux orientated drive system. The vector control is implemented in Simulink as following. The control system is in cascaded manner and consists of two major loops: the speed loop and the magnetic flux loop. The torque control loop is the minor loop of the speed control one. The references of the control loop are: the

By knowing the nameplate motor data (Fig.1), the tuning of the PI speed and flux controllers parameters can be done (see the developed Matlab file shown in the Appendix 8). The conventional vector control supposes the use of the two independent control axes: *d*-axis, along with the magnetizing flux direction, and *q* axis, in quadrature with *d*-axis. The most used field orientated control is that of the rotor orientated magnetizing flux (Fig.8). The magnetization

1

2

imr [A]

3

optimal control isq=f(t)

0

0

**4.2. Conventional vector control**

desired angular velocity for the speed loop, *ω*\*

\*

**Tuning of the speed and flux controllers**

flux and the torque can be controlled independently.

2

4

load torque [Nm]

*q*(*t*) (left to right)

for the torque loop, *i*

\* *ds*.

flux loop, *i*

6

2

4

isq [A]

354 MATLAB Applications for the Practical Engineer

6

**Figure 8.** Simulink implementation of the rotor field orientated vector control of the three-phase induction motor

**Figure 9.** Simulink implementation of the speed and flux Proportional-Integral (PI) regulators

The inputs of the speed and flux regulators Simulink block (Fig.9) are the speed reference and the flux reference; the outputs are the *q*-axis reference current, IQ\*, and *d*-axis reference current, ID\*(Fig.10). More details can be obtained by using Fig. 10. The inputs of the speed and flux regulators Simulink block (Fig.9) are speed reference and the flux reference; the outputs are: the q-axis reference current, IQ\*, and d-axis reference

current, ID\*(Fig.10). More details can be obtained by using Fig. 10.

according to Matlab file shown in the Appendix 8. **Figure 10.** PI Speed regulator of type1 (kp\_om+ki\_om/s) and Flux regulator of type 3 (kp\_df+ki\_df/s), according to Matlab file shown in the Appendix 8

Figure 10.PI Speed regulator of type1 (kp\_om+ki\_om/s) and Flux regulator of type 3 (kp\_df+ki\_df/s),

The parameters of the speed and flux regulators are delivered by using the Matlab file

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

(kp\_iqm, ki\_iqm)

(kp\_idm, ki\_idm)

w\*sigma\*Ls

1 Ts.s+1 Sample&Hold

om\_L

1 Ts.s+1

4 ID

axes, and the voltage compensation Simulink block of Fig.7

the independent control of flux and torque is obtained (Fig.12).

Modulation (PWM) transfer function block, shown in the Fig.13.

[-sqrt(2/3)\*Un, sqrt(2/3)\*Un]

Saturation

2 VD\*

1 VQ\*

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

357

3 E\_Ds

Transfer Fcn1 Saturation1

P+I/s

PI controller1

0

E\_Ds1

**Figure 12.** Detailed program of the PI control (torque (IQ\*) and flux (ID\*) control loops), the decoupling of the *d*, *q*

In the Fig.12. the detailed decoupling and the voltage compensation between *d*, *q* axes are depicted. In order to control the torque current component (IQ) and the flux current component (ID), the PI controllers are involved. By introducing the compensating voltage components,

The reference voltage components pair (VD\*, VQ\*) is used as input of the Pulse Width

Product3

Product2

P+I/s

PI controller

developed in the Appendix 8.

IQ\*

1 IQ\*

2 IQ

3 ID\*

1/2/pi freq

Gain1

By using the field orientated mechanism implemented in Simulink (Fig.11), based on the input signals (IQ\*, flux and om), the frequency of the stator voltage is delivered. By using the field orientated mechanism implemented in Simulink (Fig.11), based on the input signals (IQ\*, *flux* and *om*), the frequency of the stator voltage is delivered.

**Figure 11.** The Field Orientated Control (FOC) mechanism implemented in Simulink

stator frequency determination

developed in the Appendix 8.

1 iq\*

The parameters of the speed and flux regulators are delivered by using the Matlab file developed in the Appendix 8.

(kp\_iqm, ki\_iqm)

The inputs of the speed and flux regulators Simulink block (Fig.9) are the speed reference and the flux reference; the outputs are the *q*-axis reference current, IQ\*, and *d*-axis reference current,

current, ID\*(Fig.10). More details can be obtained by using Fig. 10.

1 Out\_1 Sum2

The inputs of the speed and flux regulators Simulink block (Fig.9) are speed reference and the flux reference; the outputs are: the q-axis reference current, IQ\*, and d-axis reference

> kp\_df Proporti onal ki\_df s Integral

Figure 10.PI Speed regulator of type1 (kp\_om+ki\_om/s) and Flux regulator of type 3 (kp\_df+ki\_df/s),

1 Out\_1 Sum2

By using the field orientated mechanism implemented in Simulink (Fig.11), based on the

1/2/pi freq

Gain1

cos(u(1)) Fcn1

The parameters of the speed and flux regulators are delivered by using the Matlab file

2 SIN(theta)

theta2

sin Trigonometric Function

1 COS(theta)

1/2/pi freq

theta3

Gain1

cos(u(1)) Fcn1

> 2 SIN(theta)

> > theta2

1 COS(theta)

theta3

input signals (IQ\*, flux and om), the frequency of the stator voltage is delivered.

**Figure 10.** PI Speed regulator of type1 (kp\_om+ki\_om/s) and Flux regulator of type 3 (kp\_df+ki\_df/s), according to

By using the field orientated mechanism implemented in Simulink (Fig.11), based on the input

sin Trigonometric Function

z 1

Unit Delay

om\_e

1 In\_1

Figure 11.The Field Orientated Control (FOC) mechanism implemented in Simulink.

ID\*(Fig.10). More details can be obtained by using Fig. 10.

kp\_om Proportional ki\_om s Integral

356 MATLAB Applications for the Practical Engineer

1 In\_1

according to Matlab file shown in the Appendix 8.

om\_e om\_slip

**Figure 11.** The Field Orientated Control (FOC) mechanism implemented in Simulink

om\_slip

stator frequency determination

signals (IQ\*, *flux* and *om*), the frequency of the stator voltage is delivered.

z 1 Unit Delay

stator frequency determination

developed in the Appendix 8.

Product1

1 iq\*

Product1

3 om

2 flux

3 om

1 iq\*

2 flux Rr\*Lm/Lr

Matlab file shown in the Appendix 8

Gain

Rr\*Lm/Lr

Gain

**Figure 12.** Detailed program of the PI control (torque (IQ\*) and flux (ID\*) control loops), the decoupling of the *d*, *q* axes, and the voltage compensation Simulink block of Fig.7

In the Fig.12. the detailed decoupling and the voltage compensation between *d*, *q* axes are depicted. In order to control the torque current component (IQ) and the flux current component (ID), the PI controllers are involved. By introducing the compensating voltage components, the independent control of flux and torque is obtained (Fig.12).

The reference voltage components pair (VD\*, VQ\*) is used as input of the Pulse Width Modulation (PWM) transfer function block, shown in the Fig.13.

<sup>0</sup> 0.5 <sup>1</sup> 1.5 <sup>0</sup>

fr frrif

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

359

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

MA 3~

Siemens 2,2 kW

torque sensor TRS200

 LIEDTKE Magnetic particle brake FRATO 2002 watercooled

ENCODER 1000ppr

timp [s]

The performances of the speed control loop are shown in the Fig. 14, in which the step rotor speed reference (om\_ref\*) is applied and the actual rotor speed (om) follows the reference with zero steady state error. The performances of the implemented flux control loop are shown in the Fig. 15. Small flux overshoot (up to 11%) and zero steady state error are obtained (Fig.15). The IM starts after *t*=0,5s (Fig.15), when the rated flux of the motor is attained (Fig.15). In order to implement the optimal control for ac drives, a test bench based on the DS1104 platform has

> uT uS uR iLc

GALVANIC ISOLATION

*DS 1104 Controller Board*  with Matlab/Simulink Interface

iLb iLa

CP1104

Tl [Nm] q

0.1

[t\* ,n\* ] PCI

**Figure 16.** The experimental set-up test bench (Gaiceanu *et al*., 2008)

ControlDesk Command Interface

**Figure 15.** The rotor flux control loop

been developed (Fig.16).

3x380V; 50Hz

0.2

0.3

flux rotor [Wb]

0.4

0.5

0.6

0.7

**Figure 13.** The Simulink models of the PWM modulator and of the equivalent *d*, *q* model of the IM

**Figure 14.** The speed control loop

**Figure 15.** The rotor flux control loop

2 ID

wr wr\*

1 IQ

w\*sigma\*Ls

1 Lssigma.s+Rs Transfer Fcn

om\_L1

1 Lssigma.s+Rs Transfer Fcn1

Product5

Product4

Kpwm Tpwm.s+1 Transfer Fcn3

angular velocity [rad/s]

**Figure 14.** The speed control loop

2 VD\*

1 VQ\*

358 MATLAB Applications for the Practical Engineer

3 E\_Ds

**Figure 13.** The Simulink models of the PWM modulator and of the equivalent *d*, *q* model of the IM

<sup>0</sup> 0.5 <sup>1</sup> 1.5 -20

time[s]

Kpwm Tpwm.s+1 Transfer Fcn2

> The performances of the speed control loop are shown in the Fig. 14, in which the step rotor speed reference (om\_ref\*) is applied and the actual rotor speed (om) follows the reference with zero steady state error. The performances of the implemented flux control loop are shown in the Fig. 15. Small flux overshoot (up to 11%) and zero steady state error are obtained (Fig.15). The IM starts after *t*=0,5s (Fig.15), when the rated flux of the motor is attained (Fig.15). In order to implement the optimal control for ac drives, a test bench based on the DS1104 platform has been developed (Fig.16).

**Figure 16.** The experimental set-up test bench (Gaiceanu *et al*., 2008)

platform.

q=0].

time).

5. Conclusions

Figure 16.The experimental set-up test bench (Gaiceanu et al., 2008).

Both types of control, i.e. conventional and optimal, have been implemented on DS1104 **Figure 17.** Real time implementation of both electrical drive systems based on the DS1104 real time platform

(Fig.17) controller board. The developed ControlDesk interface has two main functions: real time control and signals acquisition. The evaluation of the obtained experimental results reveals the system efficiency improvement through an important decreasing of the windings dissipation power and input power (Fig.18) in optimal control case. Therefore, the power balance provides an increasing of the system efficiency, thus improved thermal and capacity conditions are obtained. Both types of control, i.e. conventional and optimal, have been implemented on DS1104 (Fig. 17) controller board. The developed ControlDesk interface has two main functions: real time control and signals acquisition. The evaluation of the obtained experimental results reveals the system efficiency improvement through an important decreasing of the windings dissi‐ pation power and input power (Fig.18) in optimal control case. Therefore, the power balance provides an increasing of the system efficiency, thus improved thermal and capacity conditions are obtained.

Figure 17.Real time implementation of both electrical drive systems based on the DS1104 real time

The optimal control (7) of the electrical drive based on the induction machine (1) has been numerically simulated by using Z transform and zero order hold for a starting of a 0.75 [kW], 1480 [rpm] under a rating load torque of 4.77 [Nm]. The initial conditions of the optimal control

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

=0 *q*=0]. The free final conditions are: *t*<sup>f</sup>

In order to obtain the optimal control, the nonrecursive solution has been used to avoid the main disadvantages of the recursive solution (that is, the solution can be calculated at the current time, not backward in time; the load torque can have any shape, either linear or nonlinear, in admissible limits, due to the fact that the solution is computed at each sample

A simplified solution is proposed by avoiding the constrained optimal control problem type. By using an adequate weighting matrices S, R and Q, the magnitude constraints have been solved. In this way the states and the control signals are maintained in admissible limits.

Figure 18. Conventional and optimal control: comparative evaluation of the energy.

The energy components of the

WpJ[Ws]

WCu2[Ws]

three-phase IM

**Figure 18.** Conventional and optimal control: comparative evaluation of the energy

electromagnetic energy.

W1[Ws]

0

100

200

300

400

500

600

700

W[Ws] energy

WCu1[Ws]

W[Ws]

slightly increase in optimal control case (Fig.20).

The optimal problem shown in this chapter has been formulated for a starting period. In the Fig.18, the energy evaluation for both control methods, conventional and optimal, is depicted. According to Bellman's theorem of optimality, the optimal problem can be extended also on the braking, or reversing periods. The components of the energy of the three-phase IM (Fig.20) are: W1 the input energy, W2 the output energy, WCu1 the energy expended in the stator windings, WCu2 the energy expended in the rotor windings, WFv the viscous friction dissipated energy, WpJ the accumulated energy in inertial mass, W the

WFv[Ws]

W2[Ws]

Optimal Control

Conventional Control

From energetic point of view it could be noticed that the output energy is the same in both controls (conventional, optimal control) while the input energy decreases significantly in optimal control case due to the decreasing of the stator copper losses, even the rotor losses

=0.9[s], **x**<sup>f</sup>

=[*n*<sup>f</sup>

=1480 *q*=0].

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

361

problem are: *t*o=0[s], **x**<sup>i</sup>

**5. Conclusions**

time).

=[*n*<sup>i</sup>

The optimal control (7) of the electrical drive based on the induction machine (1) has been numerically simulated by using Z transform and zero order hold for a starting of a 0.75 [kW], 1480 [rpm] under a rating load torque of 4.77 [Nm]. The initial conditions of the optimal control problem are: to=0[s], xi=[ni=0 q=0]. The free final conditions are: tf=0.9[s], xf=[nf=1480

In order to obtain the optimal control the nonrecursive solution has been used to avoid the main disadvantages of the recursive solution (that is, the solution can be calculated at the current time, not backward in time; the load torque can have any shape: either linear or nonlinear, in admissible limits due to the fact that the solution is computed at each sample The optimal control (7) of the electrical drive based on the induction machine (1) has been numerically simulated by using Z transform and zero order hold for a starting of a 0.75 [kW], 1480 [rpm] under a rating load torque of 4.77 [Nm]. The initial conditions of the optimal control problem are: *t*o=0[s], **x**<sup>i</sup> =[*n*<sup>i</sup> =0 *q*=0]. The free final conditions are: *t*<sup>f</sup> =0.9[s], **x**<sup>f</sup> =[*n*<sup>f</sup> =1480 *q*=0].

#### **5. Conclusions**

Figure 16.The experimental set-up test bench (Gaiceanu et al., 2008).

**Figure 17.** Real time implementation of both electrical drive systems based on the DS1104 real time platform

Both types of control, i.e. conventional and optimal, have been implemented on DS1104 (Fig. 17) controller board. The developed ControlDesk interface has two main functions: real time control and signals acquisition. The evaluation of the obtained experimental results reveals the system efficiency improvement through an important decreasing of the windings dissi‐ pation power and input power (Fig.18) in optimal control case. Therefore, the power balance provides an increasing of the system efficiency, thus improved thermal and capacity conditions

reveals the system efficiency improvement through an important decreasing of the windings dissipation power and input power (Fig.18) in optimal control case. Therefore, the power balance provides an increasing of the system efficiency, thus improved thermal and capacity

The optimal control (7) of the electrical drive based on the induction machine (1) has been numerically simulated by using Z transform and zero order hold for a starting of a 0.75 [kW], 1480 [rpm] under a rating load torque of 4.77 [Nm]. The initial conditions of the optimal control problem are: to=0[s], xi=[ni=0 q=0]. The free final conditions are: tf=0.9[s], xf=[nf=1480

In order to obtain the optimal control the nonrecursive solution has been used to avoid the main disadvantages of the recursive solution (that is, the solution can be calculated at the current time, not backward in time; the load torque can have any shape: either linear or nonlinear, in admissible limits due to the fact that the solution is computed at each sample

platform.

360 MATLAB Applications for the Practical Engineer

q=0].

are obtained.

time).

conditions are obtained.

5. Conclusions

In order to obtain the optimal control, the nonrecursive solution has been used to avoid the main disadvantages of the recursive solution (that is, the solution can be calculated at the current time, not backward in time; the load torque can have any shape, either linear or nonlinear, in admissible limits, due to the fact that the solution is computed at each sample time).

A simplified solution is proposed by avoiding the constrained optimal control problem type. By using an adequate weighting matrices S, R and Q, the magnitude constraints have been solved. In this way the states and the control signals are maintained in admissible limits.

The optimal problem shown in this chapter has been formulated for a starting period. In the

extended also on the braking, or reversing periods. The components of the energy of the three-phase IM (Fig.20) are: W1 the input energy, W2 the output energy, WCu1 the energy expended in the stator windings, WCu2 the energy expended in the rotor windings, WFv the viscous friction dissipated energy, WpJ the accumulated energy in inertial mass, W the

From energetic point of view it could be noticed that the output energy is the same in both controls (conventional, optimal control) while the input energy decreases significantly in optimal control case due to the decreasing of the stator copper losses, even the rotor losses

Fig.18, the energy evaluation for both control methods, conventional and optimal, is depicted. According to Bellman's theorem of optimality, the optimal problem can be **Figure 18.** Conventional and optimal control: comparative evaluation of the energy

slightly increase in optimal control case (Fig.20).

electromagnetic energy.

Figure 18. Conventional and optimal control: comparative evaluation of the energy.

A simplified solution is proposed by avoiding the constrained optimal control problem type. By using adequate weighting matrice**s S**, **R** and **Q**, the magnitude constraints have been solved. In this way the states and the control signals are maintained in admissible limits.

**Appendix 1**

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

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

363

The optimal problem shown in this chapter has been formulated for a starting period. In the Fig.18, the energy evaluation for both control methods, conventional and optimal, is depicted. According to Bellman's theorem of optimality, the optimal problem can be extended also on the braking, or reversing periods. The components of the energy of the three-phase IM (Fig. 20) are: W1 the input energy, W2 the output energy, WCu1 the energy expended in the stator windings, WCu2 the energy expended in the rotor windings, WFv the viscous friction dissipated energy, WpJ the accumulated energy in inertial mass, W the electromagnetic energy.

From the energetic point of view, it could be noticed that the output energy is the same in both controls (conventional and optimal) while the input energy decreases significantly in the optimal control case due to the decreasing of the stator copper losses, even if the rotor losses slightly increase in the optimal control case (Fig.20).

The proposed optimal control conducts to 8% electrical energy reduction, despite of using the conventional control. The conclusion is: decreasing the input energy for the optimal control system at the same time with maintaining the output energy in both systems, classical and optimal, the efficient electric drive system is obtained.

Due to using only the negative exponentials terms, the system is stable.

The optimal control can be applied in the applications with frequent dynamic regimes. The minimization of the consumed energy conducts to an increasing of the operational period of the electrical drives, or to the electrical motor overload allowance.

The optimal control could be extended to the high phase order induction machine applications, in a decoupled *d*-*q* rotor field oriented reference frame and by using adequate control of the power inverter (Gaiceanu, 2002).

Taking into consideration that this chapter is focused on developing and proving the per‐ formances of the optimal Electrical Drive Systems (EDS), an immediate effect is represented by the positive impact on the market as an alternative for electrical efficiency. The key of the proposed chapter is the integration of the developed Matlab tool by means of the high technology, i.e. the development of the systems with microprocessor. The chapter leads to drawing up specific research and development methods in the field of electrical drives. It has been underlined that by making EDS in an efficient manner, the CO2 emission reduction will be also obtained.

The achievement of this chapter is the opportunity to attract young MSc. and PhD students in R&D activities that request both theoretical and practical knowledge. The author will have new opportunities for developing international interdisciplinary research, technological transfer, sharing of experience and cooperation in the framework of the specific Research Programmes.

A simplified solution is proposed by avoiding the constrained optimal control problem type. By using adequate weighting matrice**s S**, **R** and **Q**, the magnitude constraints have been solved.

The optimal problem shown in this chapter has been formulated for a starting period. In the Fig.18, the energy evaluation for both control methods, conventional and optimal, is depicted. According to Bellman's theorem of optimality, the optimal problem can be extended also on the braking, or reversing periods. The components of the energy of the three-phase IM (Fig. 20) are: W1 the input energy, W2 the output energy, WCu1 the energy expended in the stator windings, WCu2 the energy expended in the rotor windings, WFv the viscous friction dissipated

From the energetic point of view, it could be noticed that the output energy is the same in both controls (conventional and optimal) while the input energy decreases significantly in the optimal control case due to the decreasing of the stator copper losses, even if the rotor losses

The proposed optimal control conducts to 8% electrical energy reduction, despite of using the conventional control. The conclusion is: decreasing the input energy for the optimal control system at the same time with maintaining the output energy in both systems, classical and

The optimal control can be applied in the applications with frequent dynamic regimes. The minimization of the consumed energy conducts to an increasing of the operational period of

The optimal control could be extended to the high phase order induction machine applications, in a decoupled *d*-*q* rotor field oriented reference frame and by using adequate control of the

Taking into consideration that this chapter is focused on developing and proving the per‐ formances of the optimal Electrical Drive Systems (EDS), an immediate effect is represented by the positive impact on the market as an alternative for electrical efficiency. The key of the proposed chapter is the integration of the developed Matlab tool by means of the high technology, i.e. the development of the systems with microprocessor. The chapter leads to drawing up specific research and development methods in the field of electrical drives. It has been underlined that by making EDS in an efficient manner, the CO2 emission reduction will

The achievement of this chapter is the opportunity to attract young MSc. and PhD students in R&D activities that request both theoretical and practical knowledge. The author will have new opportunities for developing international interdisciplinary research, technological transfer, sharing of experience and cooperation in the framework of the specific Research

In this way the states and the control signals are maintained in admissible limits.

energy, WpJ the accumulated energy in inertial mass, W the electromagnetic energy.

slightly increase in the optimal control case (Fig.20).

optimal, the efficient electric drive system is obtained.

power inverter (Gaiceanu, 2002).

362 MATLAB Applications for the Practical Engineer

be also obtained.

Programmes.

Due to using only the negative exponentials terms, the system is stable.

the electrical drives, or to the electrical motor overload allowance.

**Appendix 3**

% Weighting matrices S=[100 0; 0 0]; Q=[1\*J/1 0; 0 0.001 ]; R=[0.08] ; %bloc optimal

M=[-A B\*inv(R)\*B'; Q A']

[v,d]=eig(M); vecinit=v; valnit=d;

x=sort(val); if x(1,4)==d(2,2), aux0=d(2,2); d(2,2)=d(1,1); d(1,1)=aux0; aux9=v(:,2); v(:,2)=v(:,1); v(:,1)=aux9; elseif x(1,4)==d(3,3), aux1=d(3,3); d(3,3)=d(1,1); d(1,1)=aux1; aux2=v(:,3); v(:,3)=v(:,1); v(:,1)=aux2; elseif x(1,4)==d(4,4), aux3=d(4,4); d(4,4)=d(1,1); d(1,1)=aux3; aux4=v(:,4); v(:,4)=v(:,1); v(:,1)=aux4; end if x(1,3)==d(3,3), aux5=d(3,3); d(3,3)=d(2,2); d(2,2)=aux5; aux6=v(:,3); v(:,3)=v(:,2); v(:,2)=aux6; elseif x(1,3)==d(4,4), aux7=d(4,4); d(4,4)=d(2,2); d(2,2)=aux7; aux8=v(:,4); v(:,4)=v(:,2); v(:,2)=aux8; end if x(1,1)==d(4,4), aux10=d(4,4); d(4,4)=d(3,3); d(3,3)=aux10; aux11=v(:,4); v(:,4)=v(:,3); v(:,3)=aux11; end W=v ; %the eigenvectors D=d; %the eigenvalues

%the matrix of the canonical system

% in order to obtain only the negative eigenvalues

%eigenvectors and eigenvalues

val=[d(1,1) d(2,2) d(3,3) d(4,4)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The proper arrangement of the eigenvalues and eigenvectors % % Developed by Marian Gaiceanu % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% the arrangements of the eigenvalues and of the corresponding eigenvectors on the diagonal matrix D

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

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

365

**Appendix 2**

364 MATLAB Applications for the Practical Engineer

```
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% The proper arrangement of the eigenvalues and eigenvectors % 
% Developed by Marian Gaiceanu % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Weighting matrices
 S=[100 0; 
 0 0]; 
 Q=[1*J/1 0; 
 0 0.001 ]; 
 R=[0.08] ; 
%bloc optimal 
%the matrix of the canonical system 
M=[-A B*inv(R)*B'; 
 Q A'] 
%eigenvectors and eigenvalues 
[v,d]=eig(M); 
vecinit=v; 
valnit=d; 
% the arrangements of the eigenvalues and of the corresponding eigenvectors on the diagonal matrix D 
% in order to obtain only the negative eigenvalues
val=[d(1,1) d(2,2) d(3,3) d(4,4)]; 
x=sort(val); 
if x(1,4)==d(2,2), 
 aux0=d(2,2); 
 d(2,2)=d(1,1); 
 d(1,1)=aux0; 
 aux9=v(:,2); 
 v(:,2)=v(:,1); 
 v(:,1)=aux9; 
elseif x(1,4)==d(3,3), 
 aux1=d(3,3); 
 d(3,3)=d(1,1); 
 d(1,1)=aux1; 
 aux2=v(:,3); 
 v(:,3)=v(:,1); 
 v(:,1)=aux2; 
elseif x(1,4)==d(4,4), 
 aux3=d(4,4); 
 d(4,4)=d(1,1); 
 d(1,1)=aux3; 
 aux4=v(:,4); 
 v(:,4)=v(:,1); 
 v(:,1)=aux4; 
end
if x(1,3)==d(3,3), 
 aux5=d(3,3); 
 d(3,3)=d(2,2); 
 d(2,2)=aux5; 
 aux6=v(:,3); 
 v(:,3)=v(:,2); 
 v(:,2)=aux6; 
elseif x(1,3)==d(4,4), 
 aux7=d(4,4); 
 d(4,4)=d(2,2); 
 d(2,2)=aux7; 
 aux8=v(:,4); 
 v(:,4)=v(:,2); 
 v(:,2)=aux8; 
end
if x(1,1)==d(4,4), 
 aux10=d(4,4); 
 d(4,4)=d(3,3); 
 d(3,3)=aux10; 
 aux11=v(:,4); 
 v(:,4)=v(:,3); 
 v(:,3)=aux11; end
W=v ; %the eigenvectors 
D=d; %the eigenvalues
```
while (t<=tf)&(tau>0), tau=tf-t;

 V=W22-P\*W12; Iie=Ii-emlamt;

K=V\*emlamt\*F;

Yi1=-inv(R)\*B'\*(P\*xi);

Yi4=-inv(R)\*B'\*(T\*Yiw);

 Yi2=-inv(R)\*B'\*(K\*Yif); %comanda optima isq(i)=Yi1+Yi2+Yi4;

 %mathematical model % d, q supply voltage components usd(i)=Rs\*isd\*(1-sigma\*Ts\*(om+isqo\*kq));

%electrical speed

 % data refreshing xi=[nr(i);qu(i)]; Yi3=Yiw; Yif3=Yif;

 isa(i)=2\*isalfa/3; isb(i)=-isalfa/3+isbeta/sqrt(3); isc(i)=-isalfa/3-isbeta/sqrt(3); %three phase supply stator voltages usalfa=cos(qu(i))\*usd(i)-sin(qu(i))\*usq(i);

V\_alpha=usalfa;

V\_beta=usbeta; ualfas(i)=usalfa; ubetas(i)=usbeta; % Space Vector Modulation

 usa(i)=2\*usalfa/3; usb(i)=-usalfa/3+usbeta/sqrt(3); usc(i)=-usalfa/3-usbeta/sqrt(3); ualfas1(i)=usalfa; ubetas1(i)=usbeta; V\_alfar=Va-(1/2)\*Vb-(1/2)\*Vc; ualfar(i)=V\_alfar;

ubetar(i)=V\_betar;

end

q=qu(i); om=omu(i); i=i+1;

Appendix 6

 if t>0, wp=1\*Mr; else wp=0\*Mr; end wpr(i)=wp;

 %rotor (mechanical) speed nr(i)=30\*omu(i)/pi/p; %load torque variation limits

%

 emlamt=expm(mlambda\*tau); Z=emlamt\*E\*emlamt;

T=V\*(emlamt\*E\*Iie\*H1+Iie\*H2);

P=(W21+W22\*Z)\*(inv(W11+W12\*Z));

Yiw=Yi3\*exp(-t0/Ta)+wp\*(1-exp(-t0/Ta));

Yif=Yif3\*exp(-t0/Ta)+xf\*(1- exp(-t0/Ta));

isqo=isq(i-1);%therefore i=2 as initialization;

Me(i)=p\*km\*imr\*isq(i); %electromagnetic torque

% angle of the rotor magnetization current vector, imr

qu(i)=q+t0\*(omu(i)+isq(i)/(Tr\*imr));

 isalfa=cos(qu(i))\*isd-sin(qu(i))\*isq(i); isbeta=sin(qu(i))\*isd+cos(qu(i))\*isq(i);

usbeta=sin(qu(i))\*usd(i)+cos(qu(i))\*usq(i);

V\_betar=(sqrt(3)/2)\*Vb-(sqrt(3)/2)\*Vc;

 udr=cos(qu(i))\*V\_alfar-sin(qu(i))\*V\_betar; uqr=sin(qu(i))\*V\_alfar+cos(qu(i))\*V\_betar;

[Va Vb Vc]=svm\_mod(Vd,V\_alpha,V\_beta,t,Tpwm);

omu(i)=om\*exp(-t0\*Fv/J)+p\*(Me(i)-wp)\*(1- exp(-t0\*Fv/J))/Fv;

usq(i)=Rs\*(sigma\*Ts\*isq(i)/t0-isqo\*(-1+(sigma-1)\*Ts\*imr\*kq+sigma\*Ts\*isd\*kq)-om\*(sigma\*Ts\*isd+Ts\*(sigma-1)\*imr));

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

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

367

t=t+t0; tplot(i)=t; Vaplot(i)=Va; Vbplot(i)=Vb; Vcplot(i)=Vc; imru(i)=imr; isdu(i)=isd; isqi(i)=isqo;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Main Program % % Developed by Marian Gaiceanu % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#### **Appendix 4**

Appendix 6

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

**Appendix 4**

366 MATLAB Applications for the Practical Engineer

```
% Main Program % 
% Developed by Marian Gaiceanu % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
while (t<=tf)&(tau>0), 
 tau=tf-t; 
 emlamt=expm(mlambda*tau); 
 Z=emlamt*E*emlamt; 
 P=(W21+W22*Z)*(inv(W11+W12*Z)); 
 V=W22-P*W12; 
 Iie=Ii-emlamt; 
 T=V*(emlamt*E*Iie*H1+Iie*H2); 
 K=V*emlamt*F; 
 %
 Yi1=-inv(R)*B'*(P*xi); 
 Yiw=Yi3*exp(-t0/Ta)+wp*(1-exp(-t0/Ta)); 
 Yi4=-inv(R)*B'*(T*Yiw); 
 Yif=Yif3*exp(-t0/Ta)+xf*(1- exp(-t0/Ta)); 
 Yi2=-inv(R)*B'*(K*Yif); 
 %comanda optima 
 isq(i)=Yi1+Yi2+Yi4; 
 isqo=isq(i-1);%therefore i=2 as initialization;
 %mathematical model
 % d, q supply voltage components 
 usd(i)=Rs*isd*(1-sigma*Ts*(om+isqo*kq)); 
 usq(i)=Rs*(sigma*Ts*isq(i)/t0-isqo*(-1+(sigma-1)*Ts*imr*kq+sigma*Ts*isd*kq)-om*(sigma*Ts*isd+Ts*(sigma-1)*imr)); 
 Me(i)=p*km*imr*isq(i); %electromagnetic torque
 %electrical speed 
 omu(i)=om*exp(-t0*Fv/J)+p*(Me(i)-wp)*(1- exp(-t0*Fv/J))/Fv; 
 %rotor (mechanical) speed
 nr(i)=30*omu(i)/pi/p; 
 %load torque variation limits
 if t>0, 
 wp=1*Mr; 
 else 
 wp=0*Mr; 
 end
 wpr(i)=wp; 
 % angle of the rotor magnetization current vector, imr
 qu(i)=q+t0*(omu(i)+isq(i)/(Tr*imr)); 

 % data refreshing
 xi=[nr(i);qu(i)]; 
 Yi3=Yiw; Yif3=Yif; 
 isalfa=cos(qu(i))*isd-sin(qu(i))*isq(i); 
 isbeta=sin(qu(i))*isd+cos(qu(i))*isq(i); 
 isa(i)=2*isalfa/3; 
 isb(i)=-isalfa/3+isbeta/sqrt(3); 
 isc(i)=-isalfa/3-isbeta/sqrt(3); 
%three phase supply stator voltages
usalfa=cos(qu(i))*usd(i)-sin(qu(i))*usq(i); 
V_alpha=usalfa; 
usbeta=sin(qu(i))*usd(i)+cos(qu(i))*usq(i); 
V_beta=usbeta; 
ualfas(i)=usalfa; 
ubetas(i)=usbeta; 
% Space Vector Modulation
[Va Vb Vc]=svm_mod(Vd,V_alpha,V_beta,t,Tpwm); 
 usa(i)=2*usalfa/3; 
 usb(i)=-usalfa/3+usbeta/sqrt(3); 
 usc(i)=-usalfa/3-usbeta/sqrt(3); 
 ualfas1(i)=usalfa; 
 ubetas1(i)=usbeta; 
 V_alfar=Va-(1/2)*Vb-(1/2)*Vc; 
 ualfar(i)=V_alfar; 
 V_betar=(sqrt(3)/2)*Vb-(sqrt(3)/2)*Vc; 
 ubetar(i)=V_betar; 
 udr=cos(qu(i))*V_alfar-sin(qu(i))*V_betar; 
 uqr=sin(qu(i))*V_alfar+cos(qu(i))*V_betar; 
t=t+t0; tplot(i)=t; Vaplot(i)=Va; Vbplot(i)=Vb; Vcplot(i)=Vc; imru(i)=imr; isdu(i)=isd; isqi(i)=isqo; 
q=qu(i); om=omu(i); i=i+1; 
end
```
Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

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

369

### **Appendix 6**

**Appendix 6**

368 MATLAB Applications for the Practical Engineer

control

*%PWM %gain factor*

*%sampling time Ts=0.0001; %100us*

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Current loop:* d*- axis % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %sum of parasite (small) time constant* 

*%controller of type 1: 1+s\*tau1/(s\*tau2)*

*% controller of type 2: kp\_idm+ki\_idm/s*

*ki\_idm=1/tau2 %for Pn<10kW, else: %for perturbation rejection (Blasko) % controller of type 3: kp(1+1/sTim)*

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Similarly, Current loop:* q*- axis% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*

*%from the closed loop current transfer function* 

*%sum of the small constant time, speed loop*

*%time delay constant of the sample and hold device*

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Speed regulator of type1: kp\_om+ki\_om/s, Kessler % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Speed regulator of type2: (1+s\*tau3)/s\*tau4 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Tuning of the speed controller parameters kpt(1+sTit)/sTit % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*

*%%%%%%%%%%%%%%%%%%%%%%%%%% %Speed control loop % %%%%%%%%%%%%%%%%%%%%%%%%%%*

*a=(1+cos(fimm))/sin(fimm)* 

*kpt=J/(a\*Tetm); Tit=a^2\*Tetm* 

*kp\_om=kpt; ki\_om=kpt/Tit* 

*tau3=4\*Tetm; tau4=8\*Tetm^2/J* 

*Tpwm=Ts/2;* 

*Teim=Ts+Tpwm;* 

*tau1=Lssigma/Rs tau2=2\*Teim\*Kpwm/Rs* 

*kp\_idm=tau1/tau2* 

*Tim=10\*Teim %Blasko*

*kp\_iqm=kp\_idm; ki\_iqm=ki\_idm;* 

*Ttm=2\*Teim* 

*Tetm=Tdt+Ttm; %Imposed phase margin: fimm=45\*pi/180; % a factor*

*Tdt=Ts;* 

Appendix 8

control

The Matlab file for tuning of the PI speed and flux controllers parameters in conventional

*Kpwm=1; %daca nu se considera marimile de la traductoare in sistem unificat de tensiuni: 10V*

*ki\_idm=kp\_idm/Tim %assures a damping factor D=0.707 and a bandwidth of fbanda=fs/20* 

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Tuning of the PI speed and flux controllers parameters: Kessler variant % % Developed by Marian Gaiceanu % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*

The Matlab file for tuning of the PI speed and flux controllers parameters in conventional

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

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

371

#### **Appendix 8** Appendix 8

370 MATLAB Applications for the Practical Engineer

The Matlab file for tuning of the PI speed and flux controllers parameters in conventional control The Matlab file for tuning of the PI speed and flux controllers parameters in conventional control

```
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tuning of the PI speed and flux controllers parameters: Kessler variant % 
% Developed by Marian Gaiceanu %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%sampling time 
Ts=0.0001; %100us
%PWM
%gain factor
Kpwm=1; %daca nu se considera marimile de la traductoare in sistem unificat de tensiuni: 10V
Tpwm=Ts/2; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Current loop: d- axis %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%sum of parasite (small) time constant 
Teim=Ts+Tpwm; 
%controller of type 1: 1+s*tau1/(s*tau2)
tau1=Lssigma/Rs 
tau2=2*Teim*Kpwm/Rs 
% controller of type 2: kp_idm+ki_idm/s
kp_idm=tau1/tau2 
ki_idm=1/tau2 %for Pn<10kW, else:
%for perturbation rejection (Blasko)
% controller of type 3: kp(1+1/sTim)
Tim=10*Teim %Blasko
ki_idm=kp_idm/Tim %assures a damping factor D=0.707 and a bandwidth of fbanda=fs/20 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Similarly, Current loop: q- axis%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kp_iqm=kp_idm; 
ki_iqm=ki_idm; 

%%%%%%%%%%%%%%%%%%%%%%%%%%
%Speed control loop %
%%%%%%%%%%%%%%%%%%%%%%%%%%
%from the closed loop current transfer function 
Ttm=2*Teim 
%time delay constant of the sample and hold device
Tdt=Ts; 
%sum of the small constant time, speed loop
Tetm=Tdt+Ttm; 
%Imposed phase margin:
fimm=45*pi/180; 
% a factor
a=(1+cos(fimm))/sin(fimm) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Tuning of the speed controller parameters kpt(1+sTit)/sTit %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kpt=J/(a*Tetm); Tit=a^2*Tetm 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Speed regulator of type1: kp_om+ki_om/s, Kessler %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kp_om=kpt; ki_om=kpt/Tit 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Speed regulator of type2: (1+s*tau3)/s*tau4 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tau3=4*Tetm; tau4=8*Tetm^2/J 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
```

```
% Speed regulator of type3: kp_omk+ki_omk/s %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kp_omk=tau3/tau4 %Kessler
ki_omk=1/tau4 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Flux loop (d axis)_ approximated mathematical model 1/(s*Tr/Lm): Kessler %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flux regulator of type1: (1+s*tau3f)/s*tau4f
Tr=Lr/Rr;%rotor time constant
tau3f=4*Tetm 
tau4f=8*Tetm^2/(Tr/Lm) 
% Flux regulator of type2: kp_dfk+ki_dfk/s
kp_dfk=tau3/tau4 
ki_dfk=1/tau4 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flux loop (d axis) approximated mathematical model 1/(s*Tr/Lm): Blasko%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% from closed loop d axis current transfer function
Tfm=2*Teim %the same with Ttm
%delay time constant of the flux loop sample and hold device 
Tdf=Ts; 
%time delay constant of the sample and hold device
Tefm=Tdt+Tfm; 
%Reference of the phase margin:
fimf=60*pi/180; 
% af factor
af=(1+cos(fimf))/sin(fimf) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Speed regulator of type 3:kpf(1+sTif)/sTif %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kpf=(Tr/Lm)/(a*Tefm) 
Tif=a^2*Tefm 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flux regulator of type 3: kp_df+ki_df/s %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kp_df=kpf; ki_df=kpf/Tif
```
**References**

from

[1] Athans M. & P. L. Falb (2006). Optimal Control: An Introduction to the Theory and

Matlab Based Simulation Tool of the Complete Optimal Control for Variable Speed Electrical Drives

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

373

[2] COM(2006). Action Plan for Energy Efficiency: Realising the Potential, Available

[3] Gaiceanu, M. (1997). "Optimal Control for the Variable Speed Induction Motor"- the‐ sis dissertation Master Degree, Power Electronics and Advanced Control Methods of

[4] Gaiceanu, M. (2002). Ph.D. thesis Optimal Control of the Adjustable Electrical Drive Systems with Induction Motors by using Advanced Methods, Dunarea de Jos Uni‐

[5] Gaiceanu, M., Rosu E., Paduraru R. & Dache C. (2008). Optimal Control Develop‐ ment System for Electrical Drives, *The Annals of "Dunarea de Jos" University of Galati*,

[6] Gattami A. & Rantzer A. (2005). Linear Quadratic Performance Criteria for Cascade Control, Proceedings of the 44th IEEE Conference on Decision and Control, and the

[9] Jianqiang, L., Yang, L., F., Zheng, Z., Trillion., Q. (2007). *Optimal Efficiency Control of Linear Induction Motor Drive for Linear Metro*, 2nd IEEE Conference on Industrial Ap‐ plications, 2007, ICIEA 2007, 23-25 May 2007, pp. 1981-1985, ISBN 978-1-4244-0737-8,

[11] Lorenz, R. D., Yang S.-M. (1992). "Efficiency-Optimized Flux Trajectories for Closed-Cycle Operation of Field-Orientation Induction Machines Drives", IEEE Trans. Ind.

[12] Lorenz, R. D., Yang, S.M. (1998). "Efficiency-Optimized Flux Trajectories for Closed Cycle Operation of Field Oriented Induction Machine Drives", Industry Applications Society Annual Meeting, 1988., Conference Record of the 1988 IEEE, pp 457-462, 1998

[13] Matinfar, M., K.Hashtrudi-Zaad (2005). Optimization-based Robot Compliance Con‐ trol: Geometric and Linear Quadratic Approaches, International Journal of Robotics Research, ISSN 0278-3649, Sage Publications, Inc., Vol.24, pp. 645-656, August 2005

European Control Conference 2005 Seville, Spain, December 12-15, 2005

the Electromechanical Systems, Dunarea de Jos University of Galati, June 1997

Its Applications. ISBN: 0486453286, USA, Dover, 2006

FASCICLE III, Vol.31, No.1, ISSN 1221-454X, pp. 5-10.

[7] http://ec.europa.eu/energy/action\_plan\_energy\_efficiency/doc/

[10] Leonhard, W.: *Control of Electrical Drives*, Springer-Verlag, Berlin 1996

Appl.,Vol. 28, No. 3, May/Jun. 1992, pp. 574-580.

versity of Galati, April, 2002, Romania

com\_2006\_0545\_en.pdf

2007

[8] http://www.automation.siemens.com/

#### **Acknowledgements**

This work was supported by a grant of the Romanian National Authority for Scientific Research, CNDI– UEFISCDI, project number PN-II-PT-PCCA-2011-3.2-1680.

#### **Author details**

Marian Gaiceanu\*

Dunarea de Jos University of Galati, Romania

#### **References**

*% Speed regulator of type3: kp\_omk+ki\_omk/s % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*

*% Flux regulator of type1: (1+s\*tau3f)/s\*tau4f*

*% from closed loop* d *axis current transfer function*

*%time delay constant of the sample and hold device*

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Speed regulator of type 3:kpf(1+sTif)/sTif % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Flux regulator of type 3: kp\_df+ki\_df/s % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*

*%delay time constant of the flux loop sample and hold device* 

*% Flux regulator of type2: kp\_dfk+ki\_dfk/s*

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Flux loop (d axis)\_ approximated mathematical model 1/(s\*Tr/Lm): Kessler % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Flux loop (*d *axis) approximated mathematical model 1/(s\*Tr/Lm): Blasko% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*

This work was supported by a grant of the Romanian National Authority for Scientific

Research, CNDI– UEFISCDI, project number PN-II-PT-PCCA-2011-3.2-1680.

*kp\_omk=tau3/tau4 %Kessler*

*Tr=Lr/Rr;%rotor time constant*

*Tfm=2\*Teim %the same with Ttm*

*%Reference of the phase margin:*

*af=(1+cos(fimf))/sin(fimf)* 

*kpf=(Tr/Lm)/(a\*Tefm) Tif=a^2\*Tefm* 

*kp\_df=kpf; ki\_df=kpf/Tif*

Dunarea de Jos University of Galati, Romania

**Acknowledgements**

**Author details**

Marian Gaiceanu\*

*tau4f=8\*Tetm^2/(Tr/Lm)* 

*ki\_omk=1/tau4* 

372 MATLAB Applications for the Practical Engineer

*tau3f=4\*Tetm* 

*Tdf=Ts;* 

*Tefm=Tdt+Tfm;* 

*fimf=60\*pi/180; % af factor*

*kp\_dfk=tau3/tau4 ki\_dfk=1/tau4* 


[14] Mendes, E., A. Baba, A. Razek, (1995). "Losses Minimization of a Field Oriented Con‐ trolled Induction Machine", Proceed. of IEE Electrical Machines and Drives Conf., Sep. 1995, pp. 310-314.

**Chapter 13**

**Modeling of Control Systems**

Didier López-Mancilla, Francisco G. Peña-Lecona, Miguel Mora-González and Jesús Muñoz Maciel

In studying control systems the reader must be able to model dynamical systems in mathe‐ matical terms and analyze their dynamic characteristics. This section provides an introduction to basic concepts and methodologies on modeling control systems. It also introduces some fundamentals to solve realistic models used basically in mechanical, electrical, thermal, economic, biological, and so on. A mathematical model is composed by a set of equations that describe a real system accurately, or at least fairly well. However a mathematical model is not unique for a given system, and the system under study can be represented in many different ways, in the same way a mathematical model can be used to represent different systems. Algorithms used to solve the set of equations that represent a control system require a great amount of programming instructions. Matlab is a tool that simplifies and accelerates such algorithms allowing to modeling a great variety of control systems in a very elegant way [1]. There are special Matlab toolbox useful to solve different control systems, in particular Control System Toolbox (included in MATLAB Version 7.8.0.347 (R2009a)): Creating linear models, data extraction, time-domain analysis, frequency-domain analysis, state space models, etc.

Some of these are used throughout the chapter to facilitate algorithm development.

This chapter is organized as follows; the section 1 is an introduction to modeling control systems. In section 2, some applications using electrical circuits for series and parallel circuits are given. Section 3, a second order analysis is presented. In section 4, control systems for mechanical vibrations are analyzed. Optical control systems are given in section 5, given a perspective with two basic systems: laser diode, and optical fiber. Some conclusions are

> © 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.

Additional information is available at the end of the chapter

Roger Chiu, Francisco J. Casillas,

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

**1. Introduction**

presented in section 6.


### **Chapter 13**

## **Modeling of Control Systems**

Roger Chiu, Francisco J. Casillas, Didier López-Mancilla, Francisco G. Peña-Lecona, Miguel Mora-González and Jesús Muñoz Maciel

Additional information is available at the end of the chapter

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

#### **1. Introduction**

[14] Mendes, E., A. Baba, A. Razek, (1995). "Losses Minimization of a Field Oriented Con‐ trolled Induction Machine", Proceed. of IEE Electrical Machines and Drives Conf.,

[15] Rosu, E. (1985). A Solution of Optimal Problem with Quadratic Performance Criteria,

[16] Rosu, E., Gaiceanu, M., & Bivol, I. (1998). Optimal Control Strategy for AC Drives, PEMC '98, 8th International Power Electronics & Motion Control Conference, Prague,

[17] Rosu, E., Gaiceanu, M., Bivol, I. (1998). "Load Torque Estimation for AC Motors", CNAE '98, The 9th Symposium on Electrical Drives, Craiova, Romania, ISBN

[18] Su, C.T. & Chiang C.L. (2004). Optimal Position/Speed Control of Induction Motor Using Improved Genetic Algorithm and Fuzzy Phase Plane Controller, Journal of

[19] Tamimi, J.M.Kh. & Jaddu, H.M. (2006). Optimal vector control of three-phase induc‐ tion machine, International Association of Science and Technology for Development, Proceedings of the 25th IASTED International Conference on Modeling, Identifica‐ tion and Control, Lanzarote, Spain, Acta Press, CA USA, ISSN 0-88986-549-3, 2006,

[20] Veerachary M. (2002). Optimal control strategy for a current source inverter fed in‐ duction motor, Computers and Electrical Engineering, Vol.28, No.4, July 2002, pp.

[21] Vukosavic, S.N., Jones, M., Levi, E., and Varga J. (2005). Rotor flux oriented control of a symmetrical six-phase induction machine, Electric Power Systems Research 75

Control and Intelligent Systems, ISSN 1480-1752, Vol.32, 2004

Sep. 1995, pp. 310-314.

374 MATLAB Applications for the Practical Engineer

Czech Republic, pp.4.160-4.165

973-9346-68-5, pp. 221-224, 1998

pp. 92-96

255-267

(2005) 142–152

Analele Universitatii "Dunarea de Jos" Galati

In studying control systems the reader must be able to model dynamical systems in mathe‐ matical terms and analyze their dynamic characteristics. This section provides an introduction to basic concepts and methodologies on modeling control systems. It also introduces some fundamentals to solve realistic models used basically in mechanical, electrical, thermal, economic, biological, and so on. A mathematical model is composed by a set of equations that describe a real system accurately, or at least fairly well. However a mathematical model is not unique for a given system, and the system under study can be represented in many different ways, in the same way a mathematical model can be used to represent different systems. Algorithms used to solve the set of equations that represent a control system require a great amount of programming instructions. Matlab is a tool that simplifies and accelerates such algorithms allowing to modeling a great variety of control systems in a very elegant way [1]. There are special Matlab toolbox useful to solve different control systems, in particular Control System Toolbox (included in MATLAB Version 7.8.0.347 (R2009a)): Creating linear models, data extraction, time-domain analysis, frequency-domain analysis, state space models, etc. Some of these are used throughout the chapter to facilitate algorithm development.

This chapter is organized as follows; the section 1 is an introduction to modeling control systems. In section 2, some applications using electrical circuits for series and parallel circuits are given. Section 3, a second order analysis is presented. In section 4, control systems for mechanical vibrations are analyzed. Optical control systems are given in section 5, given a perspective with two basic systems: laser diode, and optical fiber. Some conclusions are presented in section 6.

© 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.

#### **2. Electrical circuits**

Undoubtedly, the most classic circuit in literature is the RLC circuit. Due to its great usefulness in the study of linear systems, the model of the RLC circuit helps to understand some of the behaviors of an electrical control system. Thus, the next subsections will address to the modeling of series and parallel RLC circuits. Some behaviors will be analyzed using some tools from Matlab.

*Vi*

and

The transfer function becomes

consider an input constant force *vi*

therefore, the output signal results

(1 - *<sup>e</sup>*- *<sup>R</sup>* <sup>2</sup>*<sup>L</sup> t* cos

**Figure 2.** Block diagram for the RLC series circuit.

*vo* (*t*)=*Vi*

Laplace domain is given by

(*s*)=*RI*(*s*) + *LsI*(*s*) +

*Vo* (*s*)= <sup>1</sup>

*Vo* (*s*) *Vi*

*Vo*

Figure 2, shows the block diagram for the RLC series circuit.

In Figure 4, the step response for the RLC series circuit is plotted.

4*L* + *R* <sup>2</sup> *C* 4*L* <sup>2</sup>

(*t*)=*Vi*

(*s*) <sup>=</sup> <sup>1</sup>

(*s*)= *Vi*

In order to compute a general solution and to analyze the behavior of this RLC series circuit,

. Then *Vi*

*<sup>C</sup> <sup>t</sup>* - *<sup>R</sup>* <sup>2</sup>

Suppose a particular case for the block diagram with *R* =2*k*Ω, *L* =1*mH* and *C* =1μF. An easy way to obtain a step response for a RLC series circuit is to draw a block diagram in Simulink, from Matlab, using a source Step and a Scope to observe the response as shown in Figure 3.

*C* 4*L* + *R* <sup>2</sup>

*<sup>C</sup> <sup>e</sup>*- *<sup>R</sup>* <sup>2</sup>*<sup>L</sup> t*

1

*Cs I*(*s*), (3)

Modeling of Control Systems http://dx.doi.org/10.5772/58236 377

*Cs I*(*s*). (4)

*LCs* <sup>2</sup> <sup>+</sup> *RCs* + 1 . (5)

*<sup>s</sup>*(*LCs* <sup>2</sup> <sup>+</sup> *RCs* + 1) , (6)

sin <sup>4</sup>*<sup>L</sup>* <sup>+</sup> *<sup>R</sup>* <sup>2</sup>

*C* 4*L* <sup>2</sup>

*<sup>C</sup> <sup>t</sup>*). (7)

(*s*)=*Vi* /*s*, therefore the output signal in the

#### **Case I: RLC series circuit**

**Figure 1.** RLC series circuit.

In Figure 1, a RLC series circuit is shown. The modeling consists in to express a number of equations such that all kinds of moves are expressed in those equations. Thus, the input equation is given by

$$
\partial\_t \upsilon\_i(t) = \text{Ri}(t) + L \xrightarrow[dt]{dl} \frac{dl(t)}{dt} + \frac{1}{C} \mathbb{I}\_0^t \dot{i}(\tau) d\tau. \tag{1}
$$

Let be *vo* (*t*) the output in capacitor C. Then, the output equation is represented by

$$
\psi\_o(t) = \frac{1}{C} f\_0^t i(\tau) d\tau. \tag{2}
$$

In order to find a transfer function, the Laplace transform for both equations is applied obtaining

$$V\_f(\mathbf{s}) = RI(\mathbf{s}) + LsI(\mathbf{s}) + \frac{1}{\mathbb{C}s}I(\mathbf{s}) \tag{3}$$

and

**2. Electrical circuits**

376 MATLAB Applications for the Practical Engineer

**Case I: RLC series circuit**

**Figure 1.** RLC series circuit.

equation is given by

Let be *vo*

obtaining

*vi*

(*t*)=*Ri*(*t*) <sup>+</sup> *<sup>L</sup> di*(*t*)

*vo* (*t*)= <sup>1</sup> *C ∫* 0 *t*

from Matlab.

Undoubtedly, the most classic circuit in literature is the RLC circuit. Due to its great usefulness in the study of linear systems, the model of the RLC circuit helps to understand some of the behaviors of an electrical control system. Thus, the next subsections will address to the modeling of series and parallel RLC circuits. Some behaviors will be analyzed using some tools

In Figure 1, a RLC series circuit is shown. The modeling consists in to express a number of equations such that all kinds of moves are expressed in those equations. Thus, the input

> *dt* + 1 *C ∫* 0 *t*

(*t*) the output in capacitor C. Then, the output equation is represented by

In order to find a transfer function, the Laplace transform for both equations is applied

*i*(*τ*)*dτ*. (1)

*i*(*τ*)*dτ*. (2)

$$V\_o(s) = \frac{1}{\mathbb{C}s} I\text{(s)}.\tag{4}$$

The transfer function becomes

$$\frac{V\_o(s)}{V\_i(s)} = \frac{1}{LCs^2 + RCs + 1} \,\, . \tag{5}$$

In order to compute a general solution and to analyze the behavior of this RLC series circuit, consider an input constant force *vi* (*t*)=*Vi* . Then *Vi* (*s*)=*Vi* /*s*, therefore the output signal in the Laplace domain is given by

$$V\_{\phi} \text{(s)} = \frac{V\_i}{s \left(LCs^2 + RCs + 1\right)} \text{ \text{\textquotedblleft}n\text{\textquotedblright}} \tag{6}$$

therefore, the output signal results

$$w\_o(t) = V\_i \left(1 - e^{-\frac{R}{2L}t} \cos\sqrt{\frac{4L + \star R^2 C}{4L^2 C}}t - \sqrt{\frac{R^2 C}{4L + \star R^2 C}}e^{-\frac{R}{2L}t}\sin\sqrt{\frac{4L + \star R^2 C}{4L^2 C}}t\right). \tag{7}$$

Figure 2, shows the block diagram for the RLC series circuit.

**Figure 2.** Block diagram for the RLC series circuit.

Suppose a particular case for the block diagram with *R* =2*k*Ω, *L* =1*mH* and *C* =1μF. An easy way to obtain a step response for a RLC series circuit is to draw a block diagram in Simulink, from Matlab, using a source Step and a Scope to observe the response as shown in Figure 3. In Figure 4, the step response for the RLC series circuit is plotted.

In Figure 5, a RLC parallel circuit is shown. The equation of current could be represented by

*vL* (*τ*)*dτ* + *C*

(*τ*)*dτ* + *C*

(*s*) + *CsVo*

To compute a general solution and to analyze the behavior of this RLC parallel circuit, consider an input force of constant current *i*(*t*)= *I*. Then *I*(*s*)= *I* /*s*, therefore the output signal in the

*d vo* (*t*)

*d vC*(*t*)

(*t*)=*vR*(*t*)=*vL* (*t*)=*vC*(*t*), (9)

*dt* . (8)

Modeling of Control Systems http://dx.doi.org/10.5772/58236 379

*dt* . (10)

(*s*). (11)

*RLCs* <sup>2</sup> <sup>+</sup> *Ls* <sup>+</sup> *<sup>R</sup>* . (12)

*RLCs* <sup>2</sup> <sup>+</sup> *Ls* <sup>+</sup> *<sup>R</sup> I*, (13)

<sup>4</sup>*<sup>R</sup>* <sup>2</sup>*<sup>L</sup> <sup>C</sup>* <sup>2</sup> *<sup>t</sup>*. (14)

*<sup>i</sup>*(*t*)= *vR*(*t*) *<sup>R</sup>* + 1 *<sup>L</sup> ∫* 0 *t*

*vo*

*<sup>i</sup>*(*t*)= *vo* (*t*) *<sup>R</sup>* + 1 *<sup>L</sup> ∫* 0 *t vo*

In order to find the transfer function, the Laplace transform is

*<sup>I</sup>*(*s*)= *Vo*

Then, the transfer function results

Laplace domain is given by

and the solution signal results

*vo*

(*s*) *<sup>R</sup>* + 1 *Ls Vo*

*Vo* (*s*) *<sup>I</sup>* (*s*) <sup>=</sup> *RLs*

*Vo*

(*t*)=2*IR <sup>L</sup>* 4*R* <sup>2</sup>

the eigenvalues, and the solution to state space equation.

(*s*)= *RL*

*<sup>C</sup>* - *<sup>L</sup> <sup>e</sup>*- <sup>1</sup> <sup>2</sup>*RC t*

Figure 6, shows the corresponding block diagram for the RLC parallel circuit.

sin <sup>4</sup>*<sup>R</sup>* <sup>2</sup>

In order to analyze the stability using Matlab, consider a RLC parallel circuit with *R* =1*k*Ω, *L* =1*mH* and *C* =1μF. The aim is to obtain the transfer function, a state space representation,

The coefficients of the transfer function are stored in two vectors *num* and *den* for the numerator and denominator, respectively. The following solution uses the instructions *tf*, *tf2ss,* from Control Systems Toolbox. Window a) shows this algorithm obtaining the transfer function

*C* - *L*

The output voltage through capacitor C, establishes the equalities

such that the mathematical model can be written in only one equation

**Figure 3.** Block diagram using Simulink for a step response.

**Figure 4.** Step response for RLC series circuit using Simulink.

**Case 2: RLC Parallel circuit**

**Figure 5.** RLC parallel circuit.

In Figure 5, a RLC parallel circuit is shown. The equation of current could be represented by

$$\dot{q}(t) = \frac{v\_{\mathbb{R}}(t)}{R} + \frac{1}{L} \zeta\_0^t \upsilon\_L \text{ (\text{\textquotedblleft}t\text{\textquotedblright})}{}\_{L} \text{ (\text{\textquotedblleft}t\text{\textquotedblright})} \text{d}\tau + \mathcal{C} \frac{d \, v\_{\mathbb{C}}(t)}{dt} . \tag{8}$$

The output voltage through capacitor C, establishes the equalities

$$
\upsilon\_o(t) = \upsilon\_R(t) = \upsilon\_L \tag{9}
$$

such that the mathematical model can be written in only one equation

$$\dot{q}(t) = \frac{v\_o(t)}{R} + \frac{1}{L} \zeta\_0^t v\_o(\tau) d\tau + C \frac{d|v\_o(t)|}{dt} \,. \tag{10}$$

In order to find the transfer function, the Laplace transform is

$$I\{\mathbf{s}\} = \frac{V\_o\{\mathbf{s}\}}{R} + \frac{1}{Ls}V\_o\{\mathbf{s}\} + \mathbf{C}sV\_o\{\mathbf{s}\}.\tag{11}$$

Then, the transfer function results

**Figure 3.** Block diagram using Simulink for a step response.

378 MATLAB Applications for the Practical Engineer

**Figure 4.** Step response for RLC series circuit using Simulink.

**Case 2: RLC Parallel circuit**

**Figure 5.** RLC parallel circuit.

$$\frac{V\_o(s)}{I(s)} = \frac{RLs}{RLCs\,\,^2 \,\star \, Ls \,\, \star R} \,\,. \tag{12}$$

To compute a general solution and to analyze the behavior of this RLC parallel circuit, consider an input force of constant current *i*(*t*)= *I*. Then *I*(*s*)= *I* /*s*, therefore the output signal in the Laplace domain is given by

$$V\_o \text{(s)} = \frac{RL}{RLCs^2 \star Ls \star R} I\_\prime \tag{13}$$

and the solution signal results

$$w\_o(t) = 2IR\sqrt{\frac{L}{4\pi^2C \cdot L}} e^{-\frac{1}{2RC}t} \sin\sqrt{\frac{4R^2C \cdot L}{4R^2L \cdot C^2}}t. \tag{14}$$

Figure 6, shows the corresponding block diagram for the RLC parallel circuit.

In order to analyze the stability using Matlab, consider a RLC parallel circuit with *R* =1*k*Ω, *L* =1*mH* and *C* =1μF. The aim is to obtain the transfer function, a state space representation, the eigenvalues, and the solution to state space equation.

The coefficients of the transfer function are stored in two vectors *num* and *den* for the numerator and denominator, respectively. The following solution uses the instructions *tf*, *tf2ss,* from Control Systems Toolbox. Window a) shows this algorithm obtaining the transfer function

To obtain the solution to the state space equation using Matlab, it is necessary to construct an object *S*, of class '*sym*', from *A*, using the expression *S=sym*(*A*), in order of using the instruction *ilaplace*, from Symbolic Math Toolbox. It is well known that the solution to state equation is

To obtain the solution to the state space equation using Matlab, it is necessary to construct an object S, of class 'sym', from A, using the expression S = sym(A), in order of using the instruction ilaplace, from Symbolic Math Toolbox. It is well known that the solution to state equation is given by

Then, window d) gives a solution using symbolic mathematics and the instruction *ilaplace*,

Then, window d) gives a solution using symbolic mathematics and the

Note that, all solutions have a negative exponential or they are written explicitly as denominators with positive exponential, which mean that the system is stable. This kind of results could be useful for the next sections for

Note that, all solutions have a negative exponential or they are written explicitly as denomi‐ nators with positive exponential, which mean that the system is stable. This kind of results could be useful for the next sections for mechanical vibrations and optical control systems.

[ (cos(500\*3999^(1/2)\*t) - (3999^(1/2)\*sin(500\*3999^(1/2)\*t))/3999)/exp(500\*t),

[ (3999^(1/2)\*sin(500\*3999^(1/2)\*t))/(1999500\*exp(500\*t)), (cos(500\*3999^(1/2)\*t) + (3999^(1/2)\*sin(500\*3999^(1/2)\*t))/3999)/exp(500\*t)]


In the Nature there are many examples that can be modeled with seconddegree equations. The general form is represented by a homogeneous linear differential equation with constant coefficients a, b and, c as shown below

� �� �� � �

� � ����� � �

*<sup>a</sup> <sup>x</sup>*(0) <sup>+</sup> *<sup>c</sup>*

*<sup>a</sup> <sup>x</sup>*(0) <sup>+</sup> *<sup>x</sup>*˙(0)

In the Nature there are many examples that can be modeled with second-degree equations. The general form is represented by a homogeneous linear differential equation with constant

��

� ��0� � �

�� � �� � 0, �16�

*dt* + *cx* =0, (16)

� ���� � 0, �17�

*<sup>a</sup> X* (*s*)=0, (17)

. (18)

mechanical vibrations and optical control systems.

Laplace transform of the second-order equation is ����� � ���0� � ���0� �

and solving for the dependent variable

Laplace transform of the second-order equation is

*s* <sup>2</sup>*X* (*s*) - *sx*(0) - *x*˙(0) +

*a d* <sup>2</sup>*x d t* <sup>2</sup> <sup>+</sup> *<sup>b</sup> dx*

> *b <sup>a</sup> sX* (*s*) - *<sup>b</sup>*

*<sup>X</sup>* (*s*)= *sx*(0) <sup>+</sup> *<sup>b</sup>*

*<sup>s</sup>* <sup>2</sup> <sup>+</sup> *<sup>b</sup> <sup>a</sup> <sup>s</sup>* <sup>+</sup> *<sup>c</sup> a*

3. Second Order Analysis

*x*(*t*)=L (*sI* - *A*)-1 *x*(0). (15)

Modeling of Control Systems http://dx.doi.org/10.5772/58236 381

given by

with *x*(0)= 1 1 *<sup>T</sup>* .

d)

f =

**3. Second order analysis**

coefficients a*, b* and, *c* as shown below

and solving for the dependent variable

>> s=sym('s');

Figure 7. brisi potpis

���� � ����� � �����0�. (15)

>> f=ilaplace(inv(s\*eye(2)-A))

instruction ilaplace, with ��0� � �1 1�.

**Figure 6.** Block diagram for the RLC parallel circuit.

*G*(*s*). Window b), uses the instruction *tf2ss* to convert from transfer function to state space, getting arrays A, B, C and D. Window c) illustrate the stability of the system, computing the eigenvalues from array A. Note that the eigenvalues from the system in state space is equiv‐ alent to compute the poles of transfer function *G*(*s*). Thus, the poles are placed at the negative place of the complex plane, which means that the system is stable.

To obtain the solution to the state space equation using Matlab, it is necessary to construct an object *S*, of class '*sym*', from *A*, using the expression *S=sym*(*A*), in order of using the instruction *ilaplace*, from Symbolic Math Toolbox. It is well known that the solution to state equation is given by Figure 7. brisi potpis To obtain the solution to the state space equation using Matlab, it is necessary to construct an object S, of class 'sym', from A, using the expression S =

sym(A), in order of using the instruction ilaplace, from Symbolic Math

$$\mathbf{x}(t) = \mathsf{L}\mathsf{L}\mathsf{f}(sI - A)^{-1}\mathsf{f}(0). \tag{15}$$

Then, window d) gives a solution using symbolic mathematics and the instruction *ilaplace*, with *x*(0)= 1 1 *<sup>T</sup>* . Then, window d) gives a solution using symbolic mathematics and the instruction ilaplace, with ��0� � �1 1�.

```
d) 
>> s=sym('s'); 
>> f=ilaplace(inv(s*eye(2)-A)) 
f = 
[ (cos(500*3999^(1/2)*t) - (3999^(1/2)*sin(500*3999^(1/2)*t))/3999)/exp(500*t), 
[ (3999^(1/2)*sin(500*3999^(1/2)*t))/(1999500*exp(500*t)), 
(cos(500*3999^(1/2)*t) + (3999^(1/2)*sin(500*3999^(1/2)*t))/3999)/exp(500*t)]
```
Note that, all solutions have a negative exponential or they are written explicitly as denominators with positive exponential, which mean that the system is stable. This kind of results could be useful for the next sections for mechanical vibrations and optical control systems. Note that, all solutions have a negative exponential or they are written explicitly as denomi‐ nators with positive exponential, which mean that the system is stable. This kind of results could be useful for the next sections for mechanical vibrations and optical control systems.

> In the Nature there are many examples that can be modeled with seconddegree equations. The general form is represented by a homogeneous linear

#### differential equation with constant coefficients a, b and, c as shown below **3. Second order analysis**

*G*(*s*). Window b), uses the instruction *tf2ss* to convert from transfer function to state space, getting arrays A, B, C and D. Window c) illustrate the stability of the system, computing the eigenvalues from array A. Note that the eigenvalues from the system in state space is equiv‐ alent to compute the poles of transfer function *G*(*s*). Thus, the poles are placed at the negative

place of the complex plane, which means that the system is stable.

**Figure 6.** Block diagram for the RLC parallel circuit.

380 MATLAB Applications for the Practical Engineer

a)

1000];

s

b) >>

A = 1.0e+003 \* -1 -1000000 0.001 0 B = 1 0

D = 0 c)

>> eig(A) ans = 1.0e+004 \* -0.0500 + 3.1619i -0.0500 - 3.1619i

>> num=[1 0];

>> den=[0.000001 0.001

>> G=tf(num,den) Transfer function:


C = 1000000 0

[A,B,C,D]=tf2ss(num,den)

Laplace transform of the second-order equation is ����� � ���0� � ���0� � � � ����� � � � ��0� � � � ���� � 0, �17� In the Nature there are many examples that can be modeled with second-degree equations. The general form is represented by a homogeneous linear differential equation with constant coefficients a*, b* and, *c* as shown below

� �� �� � �

$$a\frac{d}{dt}{}^2\mathbf{x} + b\frac{d\mathbf{x}}{dt} + c\mathbf{x} = \mathbf{0},\tag{16}$$

�� � �� � 0, �16�

Laplace transform of the second-order equation is

and solving for the dependent variable

3. Second Order Analysis

���� � ����� � �����0�. (15)

$$\text{s}^2 X \text{(s)} \text{ - sx} \text{(0)} \text{ - \dot{x}} \text{(0)} + \frac{b}{a} \text{s} \text{X} \text{(s)} \text{ - \frac{b}{a} x \text{(0)} + \frac{c}{a} \text{X} \text{(s)} \text{=} \text{0},\tag{17}$$

��

and solving for the dependent variable

$$X \text{ (s)} = \frac{sx \text{ (0)} + \frac{b}{s}x \text{ (0)} + x \text{ (0)}}{s^2 + \frac{b}{s}s + \frac{c}{s}}. \tag{18}$$

The response for a unitary step input is given by the next equation

$$X \text{ (s)} = \frac{s^{\frac{2}{s}x(0) + s\left\lceil \frac{b}{a}x(0) + \frac{c}{a}\left\lceil 0 \right\rceil} \right\rceil}{s^{\frac{2}{s} + \frac{b}{a}s + \frac{c}{a}}} \bullet \frac{1}{s} \tag{19}$$

There are four possible behaviors of equation (19) [2]: overdamped, underdamped, undamped and critically damped. These results can be predicted by analyzing the characteristic equation, as shown below:

$$\begin{array}{c} \begin{array}{l} \text{positive,} & \text{overdamped} \\ \text{negative,} & \text{underdamped} \\ \text{-} \text{-} \text{ac}, & \text{undamped} \\ \text{zero,} & \text{critically damped} \end{array} \end{array} \tag{20}$$

**Figure 7.** Analysis of second-order system.

**4. Mechanical vibrations**

Structural vibrations in most cases are undesirable because they can negatively affect the appropriate operation of mechanical systems because of the increasing stresses and energy losses by dissipation especially at resonant frequencies. For these reasons, it is very important the prediction of resonances in the design of mechanical elements to prevent structural failures as fatigue after long period of time [3,4]. The main purpose of this section is to analyze and design a simple mechanical model well known as one-degree of freedom (ODOF) system and

Modeling of Control Systems http://dx.doi.org/10.5772/58236 383

In order to understand the dynamic response, in Figure 8 is illustrated the schematic of a linear

The system consists in a mass *m* connected to a spring of stiffness *k*. The other side of the spring is connected to the support. In a damped ODOF system, the coefficient *c* is the viscous damping

calculate with Matlab the frequency response to an external sinusoidal force.

damped ODOF system subjected to a harmonic excitation input.

Where the initial conditions *x*(0) and *x*˙(0) represent the initial state and the point of conver‐ gence of the response signal, respectively. The next results uses a step signal as input signal, with the instruction *step*, from Control Systems Toolbox. Thus, the graphical representation of a second-order analysis using Matlab as shown in Figure 7, can be obtained with the next code:

```
clear all;
x0=0.8;
xp0=0.35;
%%%% b*b>4*a*c
a=1;b=5;c=1;
num1=[x0 x0*b/a xp0];den1=[1 b/a c/a];
signal1=step(num1,den1);
%%%% b*b<4*a*c
a=1;b=1;c=1;
num2=[x0 x0*b/a xp0];den2=[1 b/a c/a];
signal2=step(num2,den2);
%%%% b=0
a=1;b=0;c=1;
num3=[x0 x0*b/a xp0];den3=[1 b/a c/a];
signal3=step(num3,den3);
%%%% b*b=4*a*c
a=1;b=5;c=1;
num4=[x0 x0*b/a xp0];den4=[1 b/a c/a];
signal4=step(num4,den4);
subplot(2,2,1);plot(signal1);grid;title('b*b>4*a*c, overdamped');
subplot(2,2,2);plot(signal2);grid;title('b*b<4*a*c, underdamped');
subplot(2,2,3);plot(signal3);grid;title('b=0, undamped');
subplot(2,2,4);plot(signal4);grid;title('b*b=4*a*c, critically damped');
```
**Figure 7.** Analysis of second-order system.

The response for a unitary step input is given by the next equation

*<sup>b</sup>* <sup>2</sup> - <sup>4</sup>*ac* ={

num1=[x0 x0\*b/a xp0];den1=[1 b/a c/a];

num2=[x0 x0\*b/a xp0];den2=[1 b/a c/a];

num3=[x0 x0\*b/a xp0];den3=[1 b/a c/a];

num4=[x0 x0\*b/a xp0];den4=[1 b/a c/a];

subplot(2,2,1);plot(signal1);grid;title('b\*b>4\*a\*c, overdamped'); subplot(2,2,2);plot(signal2);grid;title('b\*b<4\*a\*c, underdamped');

subplot(2,2,4);plot(signal4);grid;title('b\*b=4\*a\*c, critically damped');

subplot(2,2,3);plot(signal3);grid;title('b=0, undamped');

signal1=step(num1,den1);

signal2=step(num2,den2);

signal3=step(num3,den3);

signal4=step(num4,den4);

as shown below:

382 MATLAB Applications for the Practical Engineer

clear all; x0=0.8; xp0=0.35; %%%% b\*b>4\*a\*c a=1;b=5;c=1;

%%%% b\*b<4\*a\*c a=1;b=1;c=1;

%%%% b\*b=4\*a\*c a=1;b=5;c=1;

%%%% b=0 a=1;b=0;c=1; *<sup>X</sup>* (*s*)= *<sup>s</sup>* <sup>2</sup>*x*(0) <sup>+</sup> *<sup>s</sup> <sup>b</sup>*

*<sup>s</sup>* <sup>2</sup> <sup>+</sup> *<sup>b</sup> <sup>a</sup> <sup>s</sup>* <sup>+</sup> *<sup>c</sup> a*

*<sup>a</sup> <sup>x</sup>*(0) <sup>+</sup> *<sup>x</sup>*˙(0)

There are four possible behaviors of equation (19) [2]: overdamped, underdamped, undamped and critically damped. These results can be predicted by analyzing the characteristic equation,

> *positive*, *overdamped negative*, *underdamped* - 4*ac*, *undamped zero*, *critically damped*

Where the initial conditions *x*(0) and *x*˙(0) represent the initial state and the point of conver‐ gence of the response signal, respectively. The next results uses a step signal as input signal, with the instruction *step*, from Control Systems Toolbox. Thus, the graphical representation of a second-order analysis using Matlab as shown in Figure 7, can be obtained with the next code:

∙ 1

*<sup>s</sup>* . (19)

(20)

#### **4. Mechanical vibrations**

Structural vibrations in most cases are undesirable because they can negatively affect the appropriate operation of mechanical systems because of the increasing stresses and energy losses by dissipation especially at resonant frequencies. For these reasons, it is very important the prediction of resonances in the design of mechanical elements to prevent structural failures as fatigue after long period of time [3,4]. The main purpose of this section is to analyze and design a simple mechanical model well known as one-degree of freedom (ODOF) system and calculate with Matlab the frequency response to an external sinusoidal force.

In order to understand the dynamic response, in Figure 8 is illustrated the schematic of a linear damped ODOF system subjected to a harmonic excitation input.

The system consists in a mass *m* connected to a spring of stiffness *k*. The other side of the spring is connected to the support. In a damped ODOF system, the coefficient *c* is the viscous damping 4. Mechanical Vibrations

0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8

−0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2

Amplitude

Amplitude

<sup>0</sup> <sup>20</sup> <sup>40</sup> <sup>60</sup> <sup>80</sup> <sup>100</sup> <sup>120</sup> <sup>140</sup> 0.35

time (sec)

b=0, undamped

0 20 40 60 80 100 120 140

time (sec)

b\*b>4\*a\*c, overdamped

an external sinusoidal force.

excitation input.

Figure 7. Analysis of second-order system.

0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8

Amplitude

Amplitude

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> 0.2

<sup>0</sup> <sup>20</sup> <sup>40</sup> <sup>60</sup> <sup>80</sup> <sup>100</sup> <sup>120</sup> <sup>140</sup> 0.35

time (sec)

time (sec)

b\*b=4\*a\*c, critically damped

b\*b<4\*a\*c, underdamped

Structural vibrations in most cases are undesirable because they can negatively affect the appropriate operation of mechanical systems because of the increasing stresses and energy losses by dissipation especially at resonant frequencies. For these reasons, it is very important the prediction of resonances in the design of mechanical elements to prevent structural failures as fatigue after long period of time [3,4]. The main purpose of this section is to analyze and design a simple mechanical model well known as one-degree of freedom (ODOF) system and calculate with Matlab the frequency response to

In order to understand the dynamic response, in Figure 8 is illustrated the

**Figure 8.** Damped ODOF system

related with the loose of energy in the vibrating system. Considering a force varying harmon‐ ically as *f(t)=A sin (ωt);* where *A* is the amplitude oscillating at an angular frequency *ω* in rad/ sec; the equation of motion along the *x* direction can be written as:

$$m\ddot{\mathbf{x}} + c\dot{\mathbf{x}} + k\mathbf{x} = f(t) \tag{21}$$

Defining *r* as the frequency ratio of the natural frequency *ω*<sup>n</sup> to the angular frequency of the

1 *mω* <sup>2</sup>

Thus, it is clear that the frequency response of a damped ODOF represented by a complex number depends on *ω* of the sinusoidal input and the damping ratio *ξ*. Hence, the magnitude

(*r* <sup>2</sup> - 1)2 + (2*ξr*)2

**•** for *r ≈ 0*; that is, when *ω << ωn*; the problem can be considered as static, since the frequency of the excitation force is very small in comparison with the natural frequency obtaining also

**•** for *r=1*; namely *ω=ωn*; the system is at resonance, the amplitude of the output increases in

**•** for *r > 1*; that is, *ω >> ωn* the amplitude of the response approaches to zero and the phase

The schematic in Figure 9 is cantilever beam with a lumped mass at its free end and is used to detect mechanical resonances; the beam is fixed perpendicular to the support of a machine that can operate at different velocities, generating angular frequencies that move harmonically

Ignoring the mass of the beam, consider the following model to calculate the amplitude and

The following Matlab code plots the magnitude ratio and the phase angle of the frequency

phase angle of the dynamic displacement for angular frequencies of *0<ω<2* rad/sec.

3*x*

response in function of the angular frequency of the sinusoidal input.

, the above equation

Modeling of Control Systems http://dx.doi.org/10.5772/58236 385

(*<sup>r</sup>* <sup>2</sup> - 1) <sup>+</sup> *<sup>j</sup>*2*ξ<sup>r</sup>* . (26)

, (27)

*<sup>r</sup>* <sup>2</sup> - <sup>1</sup> . (28)

for any value of damping.

¨ + 2*x*˙ + 3*x* =*sin*(*ωt*). (29)

sinusoidal input *ω* and multiplying numerator and denominator by 1 / *ω* <sup>2</sup>

*G*( *jω*)=

ratio and the phase angle of the frequency response can be expressed as:

By inspection of equation 26 we can make some interpretations [5]:

function of the damping and the phase shift is 90o

<sup>|</sup>*G*( *<sup>j</sup>ω*)|= <sup>1</sup> / (*m<sup>ω</sup>* 2)

<sup>∅</sup>( *<sup>j</sup>ω*)=-tan-1 <sup>2</sup>*ξ<sup>r</sup>*

gives:

and

a very small phase shift.

.

*Dynamic analysis of a cantilever*

angle to 180o

along the *x* direction.

With the initial conditions equal to zero, the steady state solution of the differential equation can be obtained applying the Laplace transform of equation 21.

$$
\cos^2 \mathbf{x} \text{(s)} + \cos \mathbf{x} \text{(s)} + k \mathbf{x} \text{(s)} = F \text{(s)} \tag{22}
$$

Reordering the terms, the second order transfer function *G(s)* is obtained as:

$$\text{LC}\{s\} = \frac{\text{x}\{s\}}{\text{F}\{s\}} = \frac{1}{m s^2 + s s + k} \tag{23}$$

Replacing the terms *ω<sup>n</sup>* <sup>=</sup> (*<sup>k</sup> <sup>m</sup>*) and *<sup>ξ</sup>* <sup>=</sup> *<sup>c</sup>* 2 (*mk* ) , that represent the natural frequency and the relative damping factor of the system respectively; then *G(s)* can be written as:

$$\text{G(s)} = \frac{\frac{1}{m}}{\text{s}^2 + 2\text{\textdegree\textdegree C} + \text{\textdegree\textdegree G}^2} \tag{24}$$

When a sinusoidal input force is applied to a linear ODOF system, the response is also sinusoidal with the same oscillation frequency but with differences in amplitude and phase. In order to know the frequency response, *s* is replaced with *jω*, where *j* = -1 is the imaginary unit

$$G(j\omega) = \frac{\frac{1}{n}}{\cdots \omega^2 + \frac{1}{j} 2\xi \omega\_n \omega + \omega\_n^2}. \tag{25}$$

Defining *r* as the frequency ratio of the natural frequency *ω*<sup>n</sup> to the angular frequency of the sinusoidal input *ω* and multiplying numerator and denominator by 1 / *ω* <sup>2</sup> , the above equation gives:

$$\mathbb{E}G\{j\omega\} = \frac{\frac{1}{m\omega^2}}{\prod\left(r^2 \cdot 1\right) + \left(j2\zeta r\right)^2}.\tag{26}$$

Thus, it is clear that the frequency response of a damped ODOF represented by a complex number depends on *ω* of the sinusoidal input and the damping ratio *ξ*. Hence, the magnitude ratio and the phase angle of the frequency response can be expressed as:

$$\mathbb{E}\left|\mathbb{G}\left(j\omega\right)\right| = \frac{1 \left/ \left(m\omega^{2}\right)}{\sqrt{\left(r^{2} \cdot 1\right)^{2} \* \left(2\zeta r\right)^{2}}},\tag{27}$$

and

related with the loose of energy in the vibrating system. Considering a force varying harmon‐ ically as *f(t)=A sin (ωt);* where *A* is the amplitude oscillating at an angular frequency *ω* in rad/

c m

Figure 7. Analysis of second-order system.

0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8

Amplitude

Amplitude

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> 0.2

<sup>0</sup> <sup>20</sup> <sup>40</sup> <sup>60</sup> <sup>80</sup> <sup>100</sup> <sup>120</sup> <sup>140</sup> 0.35

time (sec)

time (sec)

b\*b=4\*a\*c, critically damped

b\*b<4\*a\*c, underdamped

Structural vibrations in most cases are undesirable because they can negatively affect the appropriate operation of mechanical systems because of the increasing stresses and energy losses by dissipation especially at resonant frequencies. For these reasons, it is very important the prediction of resonances in the design of mechanical elements to prevent structural failures as fatigue after long period of time [3,4]. The main purpose of this section is to analyze and design a simple mechanical model well known as one-degree of freedom (ODOF) system and calculate with Matlab the frequency response to

In order to understand the dynamic response, in Figure 8 is illustrated the schematic of a linear damped ODOF system subjected to a harmonic

x f

k

4. Mechanical Vibrations

0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8

−0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2

Amplitude

Amplitude

<sup>0</sup> <sup>20</sup> <sup>40</sup> <sup>60</sup> <sup>80</sup> <sup>100</sup> <sup>120</sup> <sup>140</sup> 0.35

time (sec)

b=0, undamped

0 20 40 60 80 100 120 140

time (sec)

b\*b>4\*a\*c, overdamped

an external sinusoidal force.

excitation input.

384 MATLAB Applications for the Practical Engineer

With the initial conditions equal to zero, the steady state solution of the differential equation

*mx*¨ + *cx*˙ + *kx* = *f* (*t*) (21)

*x*(*s*) + *csx*(*s*) + *kx*(*s*)= *F* (*s*) (22)

*ms* <sup>2</sup> <sup>+</sup> *cs* <sup>+</sup> *<sup>k</sup>* (23)

, that represent the natural frequency and the

<sup>2</sup> (24)

<sup>2</sup> . (25)

sec; the equation of motion along the *x* direction can be written as:

can be obtained applying the Laplace transform of equation 21.

*ms* <sup>2</sup>

Replacing the terms *ω<sup>n</sup>* <sup>=</sup> (*<sup>k</sup>*

**Figure 8.** Damped ODOF system

unit

Reordering the terms, the second order transfer function *G(s)* is obtained as:

*<sup>F</sup>* (*s*) <sup>=</sup> <sup>1</sup>

2 (*mk* )

1 *m <sup>s</sup>* <sup>2</sup> + 2*ξωn<sup>s</sup>* <sup>+</sup> *<sup>ω</sup><sup>n</sup>*

When a sinusoidal input force is applied to a linear ODOF system, the response is also sinusoidal with the same oscillation frequency but with differences in amplitude and phase. In order to know the frequency response, *s* is replaced with *jω*, where *j* = -1 is the imaginary

> 1 *m* -*ω*<sup>2</sup> <sup>+</sup> *<sup>j</sup>*2*ξωn<sup>ω</sup>* <sup>+</sup> *<sup>ω</sup><sup>n</sup>*

*<sup>G</sup>*(*s*)= *<sup>x</sup>*(*s*)

*<sup>m</sup>*) and *<sup>ξ</sup>* <sup>=</sup> *<sup>c</sup>*

*G*(*s*)=

*G*( *jω*)=

relative damping factor of the system respectively; then *G(s)* can be written as:

$$\mathcal{D}(j\omega) = -\tan^{-1}\frac{2\zeta r}{r^2 \cdot 1}.\tag{28}$$

By inspection of equation 26 we can make some interpretations [5]:


#### *Dynamic analysis of a cantilever*

The schematic in Figure 9 is cantilever beam with a lumped mass at its free end and is used to detect mechanical resonances; the beam is fixed perpendicular to the support of a machine that can operate at different velocities, generating angular frequencies that move harmonically along the *x* direction.

Ignoring the mass of the beam, consider the following model to calculate the amplitude and phase angle of the dynamic displacement for angular frequencies of *0<ω<2* rad/sec.

$$\mathbf{3\ddot{x}} + \mathbf{2\dot{x}} + \mathbf{3x} = \sin(\omega t). \tag{29}$$

The following Matlab code plots the magnitude ratio and the phase angle of the frequency response in function of the angular frequency of the sinusoidal input.

**5. Optical control systems**

∂ *Ne* <sup>∂</sup> *<sup>t</sup>* <sup>=</sup> *<sup>I</sup>*

∂ *N ph*

RCL parallel electric model proposed in [6].

Laser diode (LD) are perhaps one of the most used optical systems, in particular for those applications where the high frequency modulation is involved. A LD can be described as a device in which an electric current is converted to photons. The time dependent relation between the input electric current and the output photons are commonly described by the following pair of equations describing the time evolution of photon and carrier densities inside

*qad* - *<sup>A</sup>*(*Ne* - *Nom*)*<sup>N</sup> ph* - *Ne*

where *N ph* is the photon density. *A* is a proportionality constant, *Nom* is the minimum electron density required to obtain a positive gain, and *β* is the fraction of the spontaneous emission that is coupled the lasing mode, *τs* is the spontaneous carrier life time, *τph* is the photon life time, and *Ne* is the carrier density, *I* is the input current, *a* is the diode area, *d* is the thickness of the active region, *q* is the electron charge. These pair of equations can be modeled using the

*<sup>τ</sup>ph* + *β*

*Ne*

<sup>∂</sup> *<sup>t</sup>* <sup>=</sup> *<sup>A</sup>*(*Ne* - *Nom*)*<sup>N</sup> ph* - *<sup>N</sup> ph*

Here the LD can be modeled using the parallel circuit whose elements are:

*<sup>R</sup>* <sup>=</sup>*Rd* ( *Ith*

*<sup>L</sup>* <sup>=</sup> *Rd <sup>τ</sup>ph* ( *<sup>I</sup>* <sup>0</sup> *Ith* ) - 1

*<sup>C</sup>* <sup>=</sup> *<sup>τ</sup><sup>s</sup>*

*Rd* <sup>=</sup> <sup>2</sup>*kT <sup>q</sup>* ( <sup>1</sup> *Id*

If we consider a sinusoidal modulation in small signal around the quiescent point *I <sup>p</sup>* in the

+ *ip*sin *ωt* the equations (30) and (31) can be linearized as follows:

*<sup>τ</sup><sup>s</sup>* , (30)

Modeling of Control Systems http://dx.doi.org/10.5772/58236 387

*<sup>τ</sup><sup>s</sup>* , (31)

*<sup>I</sup>* <sup>0</sup> ) (32)

, (33)

*Rd* , (34)

). (35)

**5.1. Laser diode**

the laser medium.

and

and

form *I* = *I <sup>p</sup>*

**Figure 9.** Cantilever beam

```
m=3; c=2; k=3;
wn=sqrt(k/m);
w=[0:0.1:2];
zeta=c/(2*sqrt(m*k));
r=wn./w;
M=abs((1./(m*w.^2))./sqrt((r.^2-1).^2+(2*zeta*r).^2));
phase=-atan2((2*zeta*r),(r.^2-1))*180/pi;
subplot(2,1,1);plot(r,M)
subplot(2,1,2);plot(r,phase)
```
**Figure 10.** Frequency response.

#### **5. Optical control systems**

#### **5.1. Laser diode**

Laser diode (LD) are perhaps one of the most used optical systems, in particular for those applications where the high frequency modulation is involved. A LD can be described as a device in which an electric current is converted to photons. The time dependent relation between the input electric current and the output photons are commonly described by the following pair of equations describing the time evolution of photon and carrier densities inside the laser medium.

$$\frac{\partial N\_{\varepsilon}}{\partial t} = \frac{I}{qad} - A \{N\_{e} - N\_{om}\} N\_{ph} - \frac{N\_{e}}{T\_{s}} \, \tag{30}$$

and

m=3; c=2; k=3; wn=sqrt(k/m); w=[0:0.1:2];

386 MATLAB Applications for the Practical Engineer

r=wn./w;

**Figure 9.** Cantilever beam

**Figure 10.** Frequency response.

zeta=c/(2\*sqrt(m\*k));

subplot(2,1,1);plot(r,M) subplot(2,1,2);plot(r,phase)

M=abs((1./(m\*w.^2))./sqrt((r.^2-1).^2+(2\*zeta\*r).^2));

phase=-atan2((2\*zeta\*r),(r.^2-1))\*180/pi;

$$\frac{\partial \mathcal{N}\_{ph}}{\partial t} = A \{ \mathcal{N}\_e \text{ - } \mathcal{N}\_{om} \} \mathcal{N}\_{ph} \text{ - } \frac{\mathcal{N}\_{ph}}{\mathcal{T}\_{ph}} + \beta \frac{\mathcal{N}\_e}{\mathcal{T}\_s} \text{ } \tag{31}$$

where *N ph* is the photon density. *A* is a proportionality constant, *Nom* is the minimum electron density required to obtain a positive gain, and *β* is the fraction of the spontaneous emission that is coupled the lasing mode, *τs* is the spontaneous carrier life time, *τph* is the photon life time, and *Ne* is the carrier density, *I* is the input current, *a* is the diode area, *d* is the thickness of the active region, *q* is the electron charge. These pair of equations can be modeled using the RCL parallel electric model proposed in [6].

Here the LD can be modeled using the parallel circuit whose elements are:

$$R = R\_d \left(\frac{I\_{th}}{I^\circ}\right) \tag{32}$$

$$L = \frac{\frac{R\_d \ r\_{ph}}{I \left[\left(\frac{l\_0}{l\_{th}}\right) \cdot 1\right]}}{\left[\left(\frac{l\_0}{l\_{th}}\right) \cdot 1\right]} \tag{33}$$

$$\mathbf{C} = \frac{\tau\_s}{R\_d} \tag{34}$$

and

$$R\_d = \frac{2kT}{q} \left(\frac{1}{T\_d}\right). \tag{35}$$

If we consider a sinusoidal modulation in small signal around the quiescent point *I <sup>p</sup>* in the form *I* = *I <sup>p</sup>* + *ip*sin *ωt* the equations (30) and (31) can be linearized as follows:

$$\mathcal{Z}\_0(\mathbf{s}) = \frac{V\_0(s)}{I\_0(s)} = \frac{RLs}{RLCs\,\,^2 \,\star \,\, Ls \,\, \star R\,\,\, '}\tag{36}$$

**Figure 12.** Photon density response for different input values.

considering the phenomena mentioned above.

in Appendix A.

the drum (Figure 13).

The code that was used to generate and plot the curves showed in Figures 11 and 12, is placed

Modeling of Control Systems http://dx.doi.org/10.5772/58236 389

The manufacturing process of optical fiber usually requires several control systems, for example, the temperature control system of the furnace for melting the preform, the control system to maintain the fiber diameter constant, the system for deposit control and ultraviolet light curing of the fiber coating system and finally the traction control and fiber winding on

The problem of modeling a furnace for drawing fiber requires the consideration of variables such as the speed of stretching of the fiber, the rate of purge gas into the furnace, the furnace temperature, the diameter of the preform, etcetera. Moreover, variables such as the flow of gas within the furnace, radiative heat transport between the preform and the interior of the furnace is actually not so simple, this problem has been analyzed for a graphite furnace [7]. However, only for demonstration purposes, we will treat the problem for the linear case and without

**5.2. Temperature control system for a furnace in a fiber optic drawing process.**

or

$$\begin{split} Z\_{0}(\mathbf{s}) = \frac{V\_{0}(\mathbf{s})}{I\_{0}(\mathbf{s})} &= \frac{R\_{d}\left(\frac{\boldsymbol{\mu}}{\boldsymbol{\tau}\_{s}} + \frac{\boldsymbol{\mu}\_{\text{out}} \cdot \mathbf{n}\_{\text{out}} \cdot \mathbf{n}\_{\text{s}}^{0}}{\mathbf{r}\_{s} \cdot \mathbf{n}}\right)}{\boldsymbol{\tau}\_{\text{s}} \cdot \mathbf{s} + \frac{\boldsymbol{\mu}}{\boldsymbol{\tau}\_{s}} \left[\boldsymbol{\mu}\_{\text{ph}}^{0} + \mathbf{1} + \frac{\boldsymbol{\tau}\_{s}}{\boldsymbol{\tau}\_{\text{ph}} \cdot \left(\mathbf{1} + \boldsymbol{\tau}\_{\text{in}} \cdot \mathbf{n}\_{\text{s}}\right)^{0}}\right] + \frac{\boldsymbol{\mu}\_{\text{ph}}^{0} + \beta \left[\boldsymbol{\mu}\_{\text{s}}^{0} \left(\mathbf{1} + \frac{\boldsymbol{\tau}\_{\text{out}}}{\boldsymbol{\tau}\_{\text{ph}} \cdot \mathbf{n}}\right) \cdot \mathbf{n}\_{\text{s}}\right]}{\left(\beta \mathbf{1}\right) \left(\mathbf{1} + \boldsymbol{\tau}\_{\text{out}} \cdot \mathbf{n}\_{\text{s}}\right)^{0}}. \tag{37}$$

Figure 11, shows a comparative curve for a sinusoidal input and the photon density output, here; we observe a phase difference between the input and output.

**Figure 11.** Comparative curve for sinusoidal Input current *vs* photon density.

Figure 12, shows the photon density response for different values of input signal.

**Figure 12.** Photon density response for different input values.

*Z*0

(*s*) <sup>=</sup> *Rd*

here; we observe a phase difference between the input and output.

**Figure 11.** Comparative curve for sinusoidal Input current *vs* photon density.

Figure 12, shows the photon density response for different values of input signal.


or

*Z*0 (*s*)= *<sup>V</sup>*<sup>0</sup> (*s*) *I*0

388 MATLAB Applications for the Practical Engineer

(*s*)= *<sup>V</sup>*0(*s*)

*<sup>I</sup>*0(*s*) <sup>=</sup> *RLs*

( *is <sup>τ</sup><sup>s</sup>* +

*τs <sup>τ</sup> ph* (1 <sup>+</sup> *nom* - *ne*

Figure 11, shows a comparative curve for a sinusoidal input and the photon density output,

<sup>1</sup> <sup>+</sup> *nom* - *ne* 0 *τsτ ph* )

> 0) + *nph* <sup>0</sup> <sup>+</sup> *<sup>β</sup> ne* 0(1 + 1 *nph* <sup>0</sup> ) - *nom*

*RLCs* <sup>2</sup> <sup>+</sup> *Ls* <sup>+</sup> *<sup>R</sup>* , (36)

. (37)

*τsτ ph*

The code that was used to generate and plot the curves showed in Figures 11 and 12, is placed in Appendix A.

#### **5.2. Temperature control system for a furnace in a fiber optic drawing process.**

The manufacturing process of optical fiber usually requires several control systems, for example, the temperature control system of the furnace for melting the preform, the control system to maintain the fiber diameter constant, the system for deposit control and ultraviolet light curing of the fiber coating system and finally the traction control and fiber winding on the drum (Figure 13).

The problem of modeling a furnace for drawing fiber requires the consideration of variables such as the speed of stretching of the fiber, the rate of purge gas into the furnace, the furnace temperature, the diameter of the preform, etcetera. Moreover, variables such as the flow of gas within the furnace, radiative heat transport between the preform and the interior of the furnace is actually not so simple, this problem has been analyzed for a graphite furnace [7]. However, only for demonstration purposes, we will treat the problem for the linear case and without considering the phenomena mentioned above.

∆*QH* =*QC* - *QF* = *KT*

*KT* is the thermal capacity and *Ti*

calculated for the silicon preform is

Simplifying the heat flux in the furnace to

graphite with a thickness of 5mm is calculated with

*Req* =5× <sup>10</sup>-3 *<sup>m</sup>* (0.476 *<sup>W</sup>*

Then, applying the Laplace transform to Equation 38 we obtain

and 300 mm long.

in Figure 14

in Figure 15.

Where *QC* is the energy provided by the controller, *QF* is the gas heat flow into the furnace,

We estimate the thermal capacity *KT* for a preform of 100 mm of diameter and 450 mm long. We propose the active cavity of the furnace as a cylindrical shape with 120 mm in diameter

Then, knowing the mass of the perform *MP* and silicon specific heat *CP*, the heat capacity

*KT* <sup>=</sup>*MP* <sup>×</sup>*Cp* =8.2313 *Kg* ×2.575 *Joule*

*QF* <sup>=</sup> *Ti*-*To*

The equivalent thermal resistance of the internal cylindrical wall of the furnace made of

*<sup>m</sup>* <sup>∙</sup> ° *<sup>C</sup>* <sup>∙</sup> 0.1131 *<sup>m</sup>* <sup>2</sup> ) =92.87×10-3 ° *<sup>C</sup>*

Where *Req* is the equivalent thermal resistance of the inner walls of the furnace.

*Ti* <sup>=</sup> *QC*-*QF*

Using as reference the Thermo demo of Simulink, the furnace temperature control is sketched

To run the program Fiber\_Furnace\_System.mdl is necessary to preload the characteristics of the furnace in the file Data\_furnace.m (Apendix A). The graph of the system behavior is show

is the internal temperature.

*dT <sup>i</sup>*

*Kg* <sup>∙</sup> ° *<sup>C</sup>* =21.19 *joule*

*Req* . (40)

*KT* \**<sup>s</sup>* . (42)

*dt* (38)

Modeling of Control Systems http://dx.doi.org/10.5772/58236 391

° *<sup>C</sup>* . (39)

*<sup>W</sup>* . (41)

**Figure 13.** Fiber drawing process system

In the next section, we will demonstrate the use of Simulink of Matlab for modeling a simple temperature control system for a furnace in a fiber optic drawing process.

#### **5.3. Mathematical model**

Uniform heating of the preform is due in large part to effective radiative transport inside the furnace and the control of the changes in the temperature distribution. The above with the purpose to achieve fiber softening and stretching at high speed.

Making a simple energy balance of the furnace, we have the variation of energy per unit of time (∆*QH* ) given by

$$
\Delta Q\_H = Q\_C \text{ - } Q\_F = K\_T \frac{dT\_i}{dt} \tag{38}
$$

Where *QC* is the energy provided by the controller, *QF* is the gas heat flow into the furnace, *KT* is the thermal capacity and *Ti* is the internal temperature.

We estimate the thermal capacity *KT* for a preform of 100 mm of diameter and 450 mm long. We propose the active cavity of the furnace as a cylindrical shape with 120 mm in diameter and 300 mm long.

Then, knowing the mass of the perform *MP* and silicon specific heat *CP*, the heat capacity calculated for the silicon preform is

$$K\_T = M\_P \times \text{C} \\ p = 8.2313 \text{[Kg]} \times 2.575 \left[ \frac{joule}{Kg \bullet \text{"C}} \right] = 21.19 \left[ \frac{joule}{\text{"C}} \right] \text{.} \tag{39}$$

Simplifying the heat flux in the furnace to

In the next section, we will demonstrate the use of Simulink of Matlab for modeling a simple

Uniform heating of the preform is due in large part to effective radiative transport inside the furnace and the control of the changes in the temperature distribution. The above with the

Making a simple energy balance of the furnace, we have the variation of energy per unit of

temperature control system for a furnace in a fiber optic drawing process.

purpose to achieve fiber softening and stretching at high speed.

**5.3. Mathematical model**

**Figure 13.** Fiber drawing process system

390 MATLAB Applications for the Practical Engineer

time (∆*QH* ) given by

$$Q\_F = \frac{T\_{i\text{-}} T\_o}{R\_{eq}}.\tag{40}$$

Where *Req* is the equivalent thermal resistance of the inner walls of the furnace.

The equivalent thermal resistance of the internal cylindrical wall of the furnace made of graphite with a thickness of 5mm is calculated with

$$R\_{eq} = \dots \times \frac{10^3 \text{\AA}\_m \text{\AA}}{\left[ 0.476 \left\{ \frac{W}{\text{m} \cdot ^\circ \text{C}} \right\} \cdot 0.113 \text{\AA}\_m \text{\AA} \right]} = 92.87 \times 10^{-3} \left[ \text{ $^\circ \text{C}$ } \right]. \tag{41}$$

Then, applying the Laplace transform to Equation 38 we obtain

$$T\_j = \frac{Q\_C Q\_F}{K\_I \uparrow\_s^\*} \,. \tag{42}$$

Using as reference the Thermo demo of Simulink, the furnace temperature control is sketched in Figure 14

To run the program Fiber\_Furnace\_System.mdl is necessary to preload the characteristics of the furnace in the file Data\_furnace.m (Apendix A). The graph of the system behavior is show in Figure 15.

**6. Conclusions**

**Appendix A**

12.

study and analysis of control systems.

In order to study and analyze physical systems, along of this chapter, several systems were treated, and different Matlab codes were proposed to solve, analyze and, simulate those systems. Matlab is a powerful tool, and can help to students, engineers and, scientific to develop very simple, elegant and, powerful algorithms, that allows achieving successfully the

Modeling of Control Systems http://dx.doi.org/10.5772/58236 393

The following Matlab code was used to generate and plot the curves showed in Figures 11 and

**Figure 14.** Simulink model of the perform furnace.

**Figure 15.** Graph of the system behavior. Upper: Temperature around 1500oC; Lower Controller activation for each cycle.

#### **6. Conclusions**

In order to study and analyze physical systems, along of this chapter, several systems were treated, and different Matlab codes were proposed to solve, analyze and, simulate those systems. Matlab is a powerful tool, and can help to students, engineers and, scientific to develop very simple, elegant and, powerful algorithms, that allows achieving successfully the study and analysis of control systems.

#### **Appendix A**

**Figure 15.** Graph of the system behavior. Upper: Temperature around 1500oC; Lower Controller activation for each

cycle.

**Figure 14.** Simulink model of the perform furnace.

392 MATLAB Applications for the Practical Engineer

The following Matlab code was used to generate and plot the curves showed in Figures 11 and 12.

% Determine total internal glass mass = M

% Density of preforma silice glass = 2329 kg/m^3

The authors wish to express their gratitude for financial support of this project to Departa‐ mento de Ciencias Exactas y Tecnología, of Centro Universitario de los Lagos, Universidad de

Modeling of Control Systems http://dx.doi.org/10.5772/58236 395

, Francisco J. Casillas, Didier López-Mancilla, Francisco G. Peña-Lecona,

Department of Exact Sciences and Tecnology, Universitary Center of Los Lagos, University

[2] Nise, N.S. Control Systems Engineering (3rd Ed.). CECSA, ISBN 970-24-0254-9, D.F.,

[5] Robert F. Steidel, Jr. An introduction to mechanical vibrations. JohnWiley& Sons, Inc;

[6] Katz J, Margalit S, Herder C, Wilt D, Yariv A. The Intrinsic Electrical Equivalent Cir‐ cuit of a Laser Diode. IEEE Journal of Quantum Electronics 1981; QE-17(1) 4-7

[7] Zhilong Yin and Y. Jaluria. Thermal Transport and Flow in High-Speed Optical Fiber

[3] William J. Palm III. Mechanical Vibration. John Wiley & Sons, Inc; 2007.

[4] William Bolton. Control Engineering. Prentice Hall; 1998.

Drawing; *J. Heat Transfer;* 1998; 120(4) 916-930.

% -------------------------------

Mg = pi\*inrFurnace\*inhFurnace;

Miguel Mora-González and Jesús Muñoz Maciel

of Guadalajara, Lagos de Moreno, Jalisco, México

\*Address all correspondence to: rchiu@culagos.udg.mx

[1] http://www.mathworks.com/products/matlab/

densglass = 2329;

**Acknowledgements**

Guadalajara.

Roger Chiu\*

**References**

Mexico; 2002.

1989

**Author details**

Program Data\_furnace.m. The results are showed in Fig. 15.

```
% DATA_furnace
% This script runs in conjunction with the%"Fiber_Furnace_System"
% Loading images
A= imread('Preforma_Furnace.bmp');
B= imread('Controller.bmp');
% -------------------------------
% Define the furnace geometry
% -------------------------------
% internal furnace height = 0.3 m
inhFurnace = 0.30;
% internal furnace radio = 0.16 m
inrFurnace = 0.16;
% internal wall area
inwallArea = 2*pi*inrFurnace*inhFurnace;
% -------------------------------
% Determine the equivalent thermal
% resistance for graphite wall
% -------------------------------
kg = 0.476;
LWall = .005;
Req = LWall/(kg*inwallArea);
% c = cp of silice = 2.575 J/Kg-C
c = 2.575;
% -------------------------------
% Temperature Set Controller = 1500 deg C = 1773 K
% -------------------------------
TSC = 1773;
% ha = enthalpy of air at 40 deg C = 311400 J/kg
ha = TSC*c;
% Air flow rate Mdot = 0.1 kg/sec
Mdot = 0.1;
% -------------------------------
```

```
% Determine total internal glass mass = M
% -------------------------------
% Density of preforma silice glass = 2329 kg/m^3
densglass = 2329;
Mg = pi*inrFurnace*inhFurnace;
```
#### **Acknowledgements**

The authors wish to express their gratitude for financial support of this project to Departa‐ mento de Ciencias Exactas y Tecnología, of Centro Universitario de los Lagos, Universidad de Guadalajara.

### **Author details**

Program Data\_furnace.m. The results are showed in Fig. 15.

A= imread('Preforma\_Furnace.bmp'); B= imread('Controller.bmp'); % ------------------------------- % Define the furnace geometry % ------------------------------- % internal furnace height = 0.3 m

% internal furnace radio = 0.16 m

Req = LWall/(kg\*inwallArea); % c = cp of silice = 2.575 J/Kg-C

% -------------------------------

% -------------------------------

% Air flow rate Mdot = 0.1 kg/sec

% -------------------------------

% Temperature Set Controller = 1500 deg C = 1773 K

% ha = enthalpy of air at 40 deg C = 311400 J/kg

inwallArea = 2\*pi\*inrFurnace\*inhFurnace; % ------------------------------- % Determine the equivalent thermal % resistance for graphite wall % -------------------------------

% This script runs in conjunction with the%"Fiber\_Furnace\_System"

% DATA\_furnace

394 MATLAB Applications for the Practical Engineer

% Loading images

inhFurnace = 0.30;

inrFurnace = 0.16; % internal wall area

kg = 0.476; LWall = .005;

c = 2.575;

TSC = 1773;

ha = TSC\*c;

Mdot = 0.1;

Roger Chiu\* , Francisco J. Casillas, Didier López-Mancilla, Francisco G. Peña-Lecona, Miguel Mora-González and Jesús Muñoz Maciel

\*Address all correspondence to: rchiu@culagos.udg.mx

Department of Exact Sciences and Tecnology, Universitary Center of Los Lagos, University of Guadalajara, Lagos de Moreno, Jalisco, México

#### **References**


**Chapter 14**

**Advanced Decimator Modeling with a HDL Conversion**

Today, many systems designers use software tools such as Matlab to model complex, mixedtechnology systems prior to physically building and testing the system. These tools, along with their associated toolboxes provide an effective means for the initial modeling and simulation stages in a project. Such software tools also provide the means to extract information in a relevant format to aid the physical realization [see 1-3]. In this chapter we will describe the use of Simulink and its library blocks in conjunction with the HDL Coder tool in order to produce a HDL representation of the model that is, suitable for synthesis into digital logic hardware for implementation on devices such as FPGA and ASICs. We follow the idea of one single model, starting from the system level simulation to the final system integration and hardware implementation, possibly eliminating any designer interventions on the low (RTL) level. This aim can be achieved with the right concept right at the beginning of top level modeling, based

Our presentation starts with a simple example, using an existing block from the DSP System Toolbox Library [4]. We then proceed with the redesign of the same block in order to gain more control over the implementation details that will allow fulfilling specific requirements that might be needed for system integration. The presentation is accompanied with practical Simulink and HDL Coder tips that we believe will help the reader to reproduce and possibly upgrade the presented work. The two designs are first compared for equivalence in the Simulink environment, and then again with the circuit simulation of the compiled and synthesized hardware. Finally, some proposals are given to upgrade models with more

> © 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.

**in Mind**

Drago Strle and Dušan Raič

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

**1. Introduction**

Additional information is available at the end of the chapter

on the knowledge of the tool background and its limitations.

advanced features targeting the high-precision systems.
