- 1 -

DISEÑO E IMPLEMENTACIÓN DE UN SISTEMA DE CONTROL DE VELOCIDAD PARA UN MOTOR C.C.

EN ESPACIO DE ESTADOS

Martha Alicia Carranco García – A00780293 martha.carranco@emerson.com Enrique Cuan Urquizo – A00940016 enrique.cuan@hotmail.com Ruben Marbán Romero – A00807721 ruben.marban@gmail.com

Luis Ismael Minchala Ávila – A00808366 ismael.minchala@ask.ec

Resumen

El presente documento contiene el reporte teórico y el análisis de resultados de la implementación del sistema de control de velocidad de un motor de corriente continua en espacio de estados, aplicando diferentes técnicas de diseño: Realimentación de estados, Control óptimo y Control de modelo predictivo (MPC). El proceso de diseño involucra implementar el algoritmo de mínimos cuadrados por lotes para la identificación del modelo del motor, con una estructura autorregresiva con variable exógena (ARX); una vez que el modelo se encuentra validado se utiliza para el diseño de los diferentes comportamientos del controlador, puesto que se proponen estructuras de control para regulación y seguimiento de referencias arbitrarias. El proyecto se realizó como parte de los requisitos de aprobación de la cátedra de Control Digital Avanzado de la maestría en automatización del Tecnológico de Monterrey. Palabras Clave— Control, Identificación de Procesos, estados, Espacio de estados. 1. INTRODUCCIÓN

El diseño de sistemas de control más eficientes, en general, protagonizan tópicos de investigación en ésta área, y de manera particular en el diseño de controladores eficientes para el manejo de servomecanismos puesto que son elementos de uso cotidiano en sencillas aplicaciones de uso personal hasta complejas soluciones a nivel industrial. Cualquier máquina de rotación continua necesita un sistema de servocontrol, cuyos principales objetivos pueden resumirse en:

− Mantener constante la velocidad de un motor entre ciertos límites, aún cuando la carga ejercida en el eje del motor varíe. Este problema se relaciona de manera sencilla con el diseño de un regulador, en el campo del control. Aplica un criterio idéntico en el caso de control de la posición de la flecha del motor.

− Variar la velocidad y la carga de acuerdo a requerimientos externos del sistema. Esto también es conocido como seguimiento de una referencia.

El presente documento está organizado de la siguiente manera: En la sección 2 se presenta la metodología

utilizada para la identificación del modelo de la planta a través del método de identificación de mínimos cuadrados. La sección 3, por otro lado detalla el diseño y la implementación del sistema de control en espacio de estados utilizando un vector de realimentación y dinámica integral adicional para lidiar contra perturbaciones. En la sección número 4 se presentan los diseños y resultados obtenidos en sistemas de

- 2 -

seguimiento de referencias arbitrarias. El diseño del controlador óptimo, al igual que su aplicación práctica se discuten en la quinta sección. En la sección 6 se describen los criterios de diseño y la implementación del controlador de modelo predictivo. La penúltima sección de éste reporte presenta la implementación de un sistema de detección de fallas a través del uso de observadores de estados y finalmente se presentan las conclusiones del trabajo. 2. IDENTIFICACIÓN DEL SISTEMA

Previo al diseño de un sistema de control, resulta vital identificar el modelo del sistema que nos interesa controlar. Identificar significa determinar el modelo de un sistema dinámico a partir de las mediciones de las entradas y salidas. La identificación de un proceso, por tanto requiere el conocimiento de la señal de estímulo y la respuesta del sistema ante ese estímulo. Se recomienda elegir como señal de excitación, una señal con un contenido rico en el espectro de frecuencias de forma que se permita cubrir todo el ancho de banda de la planta a ser identificada, siendo lo ideal aplicar una señal de ruido blanco. Sin embargo, al ser la señal de ruido blanco ideal, se propone en [1], generar una señal binaria pseudoaleatoria que se asemeja mucho en contenido frecuencial a una señal de ruido blanco. Tal señal, mencionada antes, se conoce como señal PRBS (Pseudo Randomic Binary Sequence), y será la que finalmente se aplique al proceso a identificar.

La identificación de sistemas es un enfoque experimental para determinar el modelo dinámico de un sistema, en el que se involucran 4 pasos:

− Adquisición de los datos de la entrada y salida, bajo un protocolo experimental. − Selección de la estructura del modelo. − Estimación de los parámetros del modelo. − Validación del sistema identificado

Fig. 1. Identificación del modelo dinámico del motor c.c.

En la figura 1, se ilustra la arquitectura del sistema implementado para identificar el modelo dinámico del motor de corriente continua, a través de la aplicación de una señal PRBS generada por software y que gracias a una tarjeta de adquisición de datos USB-6008 de National Instruments y un circuito de amplificación se alimenta al motor con señales de voltaje. Utilizando un canal analógico de la DAQ se realiza la medición de la respuesta de velocidad. En la figura 2 se presentan los resultados de la señal aplicada y la respuesta de velocidad. El tiempo de muestreo de la prueba fue de 50 ms.

- 3 -

Fig. 2. Respuesta del sistema ante la señal PRBS.

El método de estimación del modelo que utilizamos fue el de mínimos cuadrados por lotes, debido a su sencillez en la implementación y su confiabilidad. Se trata de ajustar un modelo ideal de estructura autorregresiva con variable exógena (ARX) a los datos obtenidos en el funcionamiento real del sistema [2].

Una vez que aplicamos el método de mínimos cuadrados a los datos reales del proceso generados a través

del experimento de aplicar la señal PRBS al sistema, validamos el modelo a través de las correlaciones cruzadas entre el error y los coeficientes de mayor o menor orden que el modelo propuesto contenga. Seguidamente se analiza la “blancura del error” a través de un periodograma y la generación de la señal de respuesta del modelo estimado en contraste con la señal de respuesta del motor, como se ilustra en la figura 3.

Fig. 3. Validación del modelo generado por el método de mínimos cuadrados.

La validación del modelo nos garantiza que éste contiene el detalle de las características dinámicas del proceso para la frecuencia de muestreo al que éste fue identificado. Finalmente, el modelo que el proceso de identificación nos entregó es el que se presenta en la ecuación 1.

( )( ) 344798.02447.1

5831.312929.262

21

+−+= −

zz

zzz

zU

zY

(1)

- 4 -

3. REALIMENTACIÓN DE ESTADOS PARA UBICACIÓN DE POLOS E N LAZO CERRADO

Un sistema es de estado completamente controlable si todas las variables de estado pueden ser controladas en un período finito, mediante alguna señal de control no restringida [3]. Por tanto, si un sistema es completamente controlable, entonces los polos del sistema en lazo cerrado pueden ubicarse en cualquier localización deseada, mediante una realimentación del estado a través de una matriz de ganancia de realimentación del estado apropiada.

La presente técnica de diseño empieza con una determinación de los polos en lazo cerrado deseados,

basados en los requisitos de respuesta transitoria y/o respuesta en frecuencia como velocidad, factor de amortiguamiento relativo o ancho de banda. La arquitectura que el sistema de control tendrá se presenta en la figura 4. De allí se puede observar que el control se basa en el conocimiento/medición de todas las variables de estado, que como se analiza en la sección 4 a través del uso de observadores, no siempre están disponibles para su medición directa. Sin embargo, la literatura formal, como [3] y [4], demuestra que el proceso de diseño de los observadores y los vectores de realimentación para ubicación de polos trabajan bien en conjunto una vez que se han sido diseñados por separado.

Fig. 4. Arquitectura del sistema de control por realimentación de estados.

Se presenta a continuación el detalle de los cálculos del vector de realimentación utilizando los polinomios

normalizados de Bessel con un tiempo de estabilización, tss = 0.5 [s], y sus respectivas pruebas en simulación y en la práctica, para cambios de referencia y perturbación.

( )( )( )( )

( )( ) ( )

( ) [ ] ( )( )

( ) motordelvelocidadkx

kx

kxky

kukx

kx

kx

kx

zz

zz

zU

zY

→∴

=

+

−=

++

+−+= −−

−−

1

2

1

2

1

2

1

21

21

01

58.31

29.26

03448.0

12447.1

1

1

3448.0245.1158.3129.26

( )( )

[ ]

−==

∆=+−=

−

06474.958.31

3032.6429.26

44591.029721.11

2

GHHM

GMJE

MC

CCTn

zzzP

(2)

- 5 -

[ ]T001887.0000269.0 −=E

( ) ( ) ( ) ( )

( ) ( )

002547.0

ˆ1

1ˆ1lim

11

1

11

1

=−=

−−−=∴

+−=+−

−

−−

→

O

zO

O

K

zzzK

krKkk

HEGG

HGIC

HxHEGx

En la figura 5 se presenta la gráfica de respuesta de simulación del sistema de control de la figura 6, utilizando los resultados obtenidos en (3) y (4). De allí se puede notar que el sistema responde acorde a las especificaciones de diseño, además de mostrar su acción correctiva ante perturbaciones que eventualmente se podrían presentar en la salida del sistema (eje del motor).

Fig. 5. Respuesta del sistema de control ante cambios de referencia y perturbación. E = [0.000269 -0.001887] 3.1. DISEÑO DEL SISTEMA DE CONTROL CON OBSERVADORES DE ESTADO

En la práctica, no todas las variables de estado están disponibles para su medición directa, incluso resulta común poder medir únicamente una variable de estado, que normalmente corresponde a la respuesta del sistema. En el caso de nuestro sistema de control, a través de un encoder incremental es posible medir la velocidad del eje del motor, lo cual corresponde a la medición directa del estado x1, definido en la ecuación (2). En casos, como estos, en los que no todos los estados se pueden medir, será necesario estimar las variables de estado que no sea posible medirlas directamente. Esta estimación se llama observación. Un observador de estados es un subsistema de control que lleva a cabo una estimación de las variables de estado, a partir de las mediciones de las variables de salida y control [3].

(3)

(4)

- 6 -

La arquitectura del sistema de control implementada para las simulaciones y el desarrollo de la aplicación práctica se presenta en la figura 6. Resulta sencillo notar en la figura 6, que el observador añade dinámica adicional al sistema, por lo que se recomienda que al elegir las ubicaciones de los polos del observador, el tiempo de estabilización sea entre 4 a 6 veces más rápido que el tiempo de estabilización utilizado en el diseño del controlador por realimentación de estados para la reubicación de polos.

Fig. 6. Arquitectura del sistema de control con el observador de estados incluido.

Se presenta a continuación el detalle de los cálculos de un vector de realimentación del observador, cuya

ubicación de polos se calculado utiliza los polinomios normalizados de Bessel con un tiempo de estabilización, tss = 0.1 [s].

( )( )

[ ]

==

∆=

+−=−

12447.1

01

01737.0102841.01

2

TO

nOC

zzzP

CGCW

JWGK

M

[ ]T327428.014186.1 −=K

En la figura 7 se presenta la gráfica de respuesta de simulación del sistema de control de la figura 7, utilizando los resultados obtenidos en (3) y (4) para el control por ubicación de polos a través del vector de realimentación y el resultado calculado en (5) para la implementación de la dinámica del observador. De allí se puede notar que el sistema responde acorde a las especificaciones de diseño, sin embargo presenta diferencias con relación a la figura 5 en el que asumimos que el sistema tenía disponible todos los estados de manera directa para su medición. Principalmente se puede notar el sobreimpulso de la respuesta de velocidad y la rápida convergencia del observador a la señal real.

(5)

- 7 -

Fig. 7. Respuesta del sistema de control ante cambios de referencia y perturbación. E = [0.000269 -0.001887]; K = [1.14186 -0.327428]

En la figura 8 se presenta la gráfica de respuesta de velocidad del motor ante cambios de referencia y perturbación, obtenida luego de haber implementado el controlador en espacio de estados utilizando LabVIEW como plataforma de desarrollo y cuya arquitectura se ilustra en la figura 6. Los resultados obtenidos en (3), (4) y (5) son considerados para la programación.

Fig. 8. Respuesta del sistema de control implementado, ante cambios de referencia y perturbación.

- 8 -

De la figura 8 ese posible concluir que el sistema responde acorde a las especificaciones, puesto que el tiempo de estabilización es ligeramente superior al diseño, siendo aproximadamente 1 segundo. Resulta sencillo notar, asimismo que el sistema es capaz de compensar las perturbaciones de carga que se presenten en el eje del motor. 4. SEGUIMIENTO A REFERENCIAS ARBITRARIAS

El problema del seguimiento a una trayectoria considera el orden del sistema para el diseño de seguimiento de una trayectoria específica. En el caso del control de posición de un motor de corriente continua, si recordamos el modelo físico, éste contienen un integrador que reduce a cero el error de estado estable ante cambios de referencia tipo escalón, sin embargo si el objetivo por ejemplo, fuese seguir una trayectoria tipo rampa, éste sistema presentará un error de estado estable diferente de cero. Se propone, por tanto incluir en el sistema de control la dinámica de la señal que se pretende seguir, generándose una arquitectura de control como la que se presenta en la figura 9.

Fig. 9. Arquitectura del sistema de control para seguimiento de referencias arbitrarias, con observador de estados.

Se presenta a continuación el detalle analítico de los cálculos para el seguimiento de trayectorias tipo:

senoidal, triangular y cuadrada. Para cada caso, se realiza el cálculo de la dinámica del observador respectivo y se presenta la simulación correspondiente. Primeramente, se realizan los cálculos del sistema de control para el seguimiento de una señal tipo rectangular. Se propone un tiempo de estabilización del sistema de control, tss = 0.5 [s], y el tiempo de estabilización del observador 10 veces menor. Asimismo, se utilizó la transformada en series de Fourier de la señal cuadrada y se utilizaron los cuatro primeros términos de la expansión para la dinámica adicional del sistema de control.

( ) ( ) ( ) ( ) ( )

( )( )( )( )( )( )( )( )( ) 19.16.207.32.307.36.29.1

7

5

3

49259

2345678

8,7

6,5

4,3

2,1

22222222

87654321

+−+−+−+−=−−−−−−−−

±=±=±=±=

++

++

++

+=

zzzzzzzzzP

ezezezezezezezez

js

js

js

js

s

d

s

c

s

b

s

as

pTpTpTpTpTpTpTpT ssssssss

ππππ

ππ

ππ

ππ

ππρ

- 9 -

[ ] 100000001

1

90211.1

61803.2

07768.3

23607.3

07768.3

61803.2

90211.1

00000001

100000090211.1

010000061803.2

001000007768.3

000100023607.3

000010007768.3

000001061803.2

000000190211.1

==

−

−

−

−

=

−

−

−

−

=

aa

aa

DC

BA

[ ]Td

ad

d

a

ad

0 0.25 1.0- 1.902 2.618- 3.077 3.236- 3.077 2.618- 1.902

0 1 0 0 0 0 0 0 0 0

0.9853- 1.985 0 0 0 0 0 0 0 0.250

0 0 0 0 0 0 0 0 0 1.0-

0 0 1 0 0 0 0 0 0 1.902

0 0 0 1 0 0 0 0 0 2.618-

0 0 0 0 1 0 0 0 0 3.077

0 0 0 0 0 1 0 0 0 3.236-

0 0 0 0 0 0 1 0 0 3.077

0 0 0 0 0 0 0 1 0 2.618-

0 0 0 0 0 0 0 0 1 1.902

=

=

=

=

B

B

BB

A

ABC

0AA

[ ][ ]0.1624- 0.1872 0 0 0 0 0 0 0 0=

=

d

d C

C

0C

Considerando un tiempo de estabilización tss = 0.5 [s]

( ) 456789101

10,9

8,76,5

4,32,1

z0005.0z002.0z0057.0z013.0z0456.0z241.0

5.0

3071.14278.105.0

9342.39657.9

5.0

6057.69898.85.0

3777.93609.7

5.0

4022.126835.4

+−+−+−=

±−=

±−=±−=

±−=±−=

zzP

jp

jp

jp

jp

jp

( )

[ ] ( )dCdddddddTT

dCCTn

T

ABABABABJE

AMJE

∆−=

∆−=−

−

19210

1

MLMMM

[ ]0246.00275.00029.00022.00011.0104101001 4-4 −×−×= −TE (6)

- 10 -

( ) ( ) ( ) ( )[ ][ ]1.13860.9790.9621.0680.9211.0680.9630.9791.1380

9210

=

⋅⋅⋅⋅=T

ssssT TTTT

ρ

ρ ρρρρ L

=

9dd

dd

d

O

AC

AC

C

WM

=

0.44 0.3 114.3 219.5 265.4 254.7 225.5 192.8 162.3 135.5

0.53 0.4 26.2 114.3 219.5 265.4 254.7 225.5 192.8 162.3

0.64 0.5 0 26.2 114.3 219.5 265.4 254.7 225.5 192.8

0.77 0.6 0 0 26.2 114.3 219.5 265.4 254.7 225.5

0.91 0.7 0 0 0 26.2 114.3 219.5 265.4 254.7

1.07 0.9 0 0 0 0 26.2 114.3 219.5 265.4

1.20 1.0 0 0 0 0 0 26.2 114.3 219.5

1.24 1.2 0 0 0 0 0 0 26.2 114.3

1.00 1.2 0 0 0 0 0 0 0 26.2

0 1.0 0 0 0 0 0 0 0 0

OW

( )( ) [ ]T

z

Oz

37.1 0 0.031- 0.051 0.051- 0.038 0.017- 0 0.011 0.009-0

0 1

=

= −

x

ρWx

( )( ) nOC

zzP

JWGK 1

102

−∆=

=

=

0.0633

0.3147

0.0031-

0.0024

0.0019-

0.0013

0.0009-

0.0005

0.0003-

0.0001

K

(7)

- 11 -

Fig. 10. Respuesta del sistema de control de seguimiento a una señal cuadrada. Respuesta de seguimiento y perturbación.

En la figura 10 se presenta el resultado de la simulación en Matlab, del sistema de control bosquejado en la figura 9, con los resultados de las ganancias (6) y (7). De allí resulta sencillo notar, que el error de estado estable es cero y que el observador converge en dos tiempos de muestreo. Se puede apreciar, asimismo la forma en que el sistema se comporta cuando una perturbación es aplicada.

Se realizan ahora los cálculos de un sistema de control para el seguimiento de una señal tipo triangular. Se propone un tiempo de estabilización del sistema de control, tss = 0.7 [s], y un tiempo de estabilización del observador de por lo menos 10 veces menor. Asimismo, se utilizó la transformada en series de Fourier de la señal triangular y se utilizaron los dos primeros términos de la expansión para la dinámica adicional del sistema de control.

( ) ( )

( )( )( )( )( ) 197.395.597.3

3

1

98

3

1

1

234

4,3

2,1

22

4321

+−+−=−−−−

±=±=

++

+=

zzzzzP

ezezezez

js

js

sss

pTpTpTpT ssss

ρ

[ ] 10001

1

97504.3

95014.5

97504.3

00001

010097504.3

001095014.5

00019705.3

==

−

−=

−

−=

aa

aa

DC

BA

- 12 -

−

−

=

=

−

−

=

=

58.31

29.26

1

97504.3

95014.5

97504.3

00000058.31

00000029.26

0000001

00010097504.3

00001095014.5

0000019705.3

d

ad

d

a

ad

B

B

BB

A

ABC

0AA

[ ][ ]0 1 0 0 0 0=

=

d

d

C

C0C

Considerando un tiempo de estabilización tss = 0.7 [s]

7.0

454.11205.77.0

4018.42613.6

7.0

53.72169.4

6,5

4,32,1

jp

jp

jp

±−=

±−=±−=

( )[ ] ( )dCddddddd

TT

dCCTn

T

ABABABABJE

AMJE

∆−=

∆−=−

−

1526

1

MLMMM

[ ]0007.00003.038.0098.025.077.0 −−−−=TE

( ) ( ) ( ) ( )[ ][ ]1622.01281.00951.00629.00313.00

5210

=

⋅⋅⋅⋅=T

ssssT TTTT

ρ

ρ ρρρρ L

=

5dd

dd

d

O

AC

AC

C

WM

( )( ) [ ]T

z

Oz

0212.000003.00011.00004.00

0 1

−−=

= −

x

ρWx

(8)

- 13 -

( )( ) nOC

zzP

JWGK 1

62

−∆=

=

=

0.1775

1.5366

0.0140-

0.0464

0.0520-

0.0197

K

Fig. 11. Respuesta del sistema de control de seguimiento a una señal triangular. Respuesta de seguimiento y perturbación.

En la figura 11 se presenta el resultado de la simulación en Matlab, del sistema de control bosquejado en la figura 9, con los resultados de las ganancias (8) y (9). De allí resulta sencillo notar, que el error de estado estable es cero y que el observador converge en seis tiempos de muestreo. Se puede apreciar, asimismo la forma en que el sistema se comporta cuando una perturbación es aplicada.

Se propone ahora el diseño de seguimiento de una señal tipo senoidal. Se utiliza un tiempo de estabilización del sistema de control, tss = 0.5 [s], y para la ubicación de las raíces del observador se considera el mismo tiempo de establecimiento que en el diseño del vector de realimentación.

( )

( )( )( ) 19754.1

1

2

2,1

22

21

+−=

−−

±=+

=

zzzP

ezez

jss

s

pTpT ss

ππ

ρ

(9)

- 14 -

[ ] 101

1

9754.1

01

19754.1

==

−=

−=

aa

aa

DC

BA

−=

=

−

−=

=

58.31

29.26

1

9754.1

03448.0058.31

12447.1029.26

0001

0019754.1

d

ad

d

a

ad

B

B

BB

A

ABC

0AA

[ ][ ]0 1 0 0=

=

d

d

C

C0C

Considerando un tiempo de estabilización tss = 0.5 [s]

( ) 0705.0549.0773.19873.26486.2

5.0

0813.39268.55.0

3142.61104.4

5.0

3142.6448.6

22345

5

4,32,1

−+−+−=

±−=

±−=±−=

zzzzzzP

jp

jp

jp

( )[ ] ( )dCddddddd

TT

dCCTn

T

ABABABABJE

AMJE

∆−=

∆−=−

−

1526

1

MLMMM

[ ]0007.00003.038.0098.025.077.0 −−−−=TE

( ) ( ) ( ) ( )[ ][ ]4540.03090.01564.00

4210

=

⋅⋅⋅⋅=T

ssssT TTTT

ρ

ρ ρρρρ L

=

3dd

dd

d

O

AC

AC

C

WM

(10)

- 15 -

( )( ) [ ]Tz

Oz

1066.000017.00019.00

0 1

−=

= −

x

ρWx

( ) ( )( ) nOC

zPzP

JWGK 1

2

−∆=

=

−

−=

1250.0

0238.1

0045.0

0048.0

K

Fig. 12. Respuesta del sistema de control de seguimiento a una señal senoidal. Respuesta de seguimiento y perturbación.

En la figura 12 se presenta el resultado de la simulación en Matlab, del sistema de control bosquejado en la figura 9, con los resultados de las ganancias (10) y (11). De allí resulta sencillo notar, que el error de estado estable es cero y que el observador converge en dos tiempos de muestreo. Se puede apreciar, asimismo la forma en que el sistema se comporta cuando una perturbación es aplicada.

Los resultados del seguimiento de referencias arbitrarias se presentan en las figuras 13 y 14. En la figura 13 se encuentra el resultado práctico de realizar el seguimiento de una señal tipo triangular, al igual que su respuesta ante una perturbación de carga en el eje del motor. Cabe resaltar que los resultados utilizados para esta implementación fueron los de las ganancias (8) y (9). La gráfica de respuesta de la figura 14, por otro lado contiene el detalle del sistema de seguimiento de una señal senoidal, para lo cual se hizo uso de los resultados (10) y (11). Resulta sencillo notar allí que el error de estado estable es cero y que el sistema es capaz de responder ante perturbaciones.

(11)

- 16 -

Fig. 13. Respuesta del sistema de control de seguimiento a una señal triangular implementado, ante cambios de referencia y perturbación.

Fig. 14. Respuesta del sistema de control de seguimiento a una señal senoidal implementado, ante cambios de referencia y perturbación.

5. DISEÑO DEL CONTROLADOR ÓPTIMO Para el diseño del controlador óptimo se utiliza la misma arquitectura del sistema de control por realimentación de estados que se que se ilustra en la figura en la figura 6. Resulta sencillo notar que el mejor desempeño en el sistema se puede conseguir a través de un vector de ganancias que optimice el funcionamiento del sistema de control. Se presentan a continuación el detalle de los cálculos de diferentes vectores óptimos, al igual que su simulación y la respuesta de velocidad del motor en la práctica utilizando el vector de mejor desempeño.

- 17 -

Vector óptimo 1

[ ]

=

−

−+=

=

=

−

−−−−

6.3 -2.9 -3.6 2.9

1.0 0 -1.00

-4430.42892.4 4430.42892.7-

-3688.32407.93689.32406.6-

1

10

01

11

H

AQA

ABBRQABBRAH

R

Q

TT

TTTT

Calculando los valores y vectores propios de H, tenemos:

0005.0

1413.0

1.7

2.2020

4

3

2

1

====

λλλλ

=

0.5918 0.4617- 0.0650 0.0005

0.3871- 0.5540 0.0789- 0.0004

0.5920 0.5400- 0.6236 0.7681-

0.3867- 0.4340 0.7750 0.6403-

P

Considerando los valores propios cuya magnitud es menor que la unidad y sus correspondientes vectores, formamos la matriz T.

∴

−−−=

−−

−−

=

I

I

Y

X

T

5918.04617.0

3871.05540.0

5920.05400.0

3867.04340.0

−−

=

−−

=

5918.04617.0

3871.05540.0

5920.05400.0

3867.04340.0

I

I

Y

X

−−

==

=

=

0003.00652.0

0002.00783.0

0005.00

01413.0

0

0

4

3

H

H

EYW

E

I

λλ

- 18 -

11 −−−= IWXBRE T

[ ]0193.00187.0=E

( ) ( ) ( )

( ) ( )014829.0

1

1ˆ1lim

1

1

1

11

1

ˆ

=

−−−=∴

+

−=+

−

−−−

→

O

zO

O

K

zzzK

krKkk

HGIC

HxHEGxG43421

Fig. 15. Respuesta del sistema de control ante cambios de referencia y perturbación, utilizando un vector óptimo, E = [0.0187 0.0193]T y observadores de estados.

Vector óptimo 2

[ ]

×=

−

−+=

=

=

−

−−−−

0.0004 -0.0003 -0.3610 0.0003

0.0001 0 -0.1000 0

-0.00440.00294.43040.0029-

-0.00370.00243.68840.0023-

101

100

10000

01

4

11

H

AQA

ABBRQABBRAH

R

Q

TT

TTTT

- 19 -

Calculando los valores y vectores propios de H, tenemos:

0

06512.0

2

44283

4

3

2

1

====

λλλλ

=

1.0000 0.6398- 0.4332 0.0625

0.0007- 0.7685 0.5400- 0.0173

0.0010 0.0011- 0.0013 0.7669-

0.0007- 0.0020 0.7216 0.6385-

P

Considerando los valores propios cuya magnitud es menor que la unidad y sus correspondientes vectores, formamos la matriz T.

∴

−−−=

−−

−−

=

I

I

Y

X

T

16398.0

0007.07685.0

0010.00011.0

0007.00020.0

−−

=

−−

=

16398.0

0007.07685.0

0010.00011.0

0007.00020.0

I

I

Y

X

−==

=

=

04166.0

005004.0

00

06512.0

0

0

4

3

H

H

EYW

E

I

λλ

11 −−−= IWXBRE T

[ ]0123.00079.0=E

( ) ( ) ( )

( ) ( )006061.0

1

1ˆ1lim

1

1

1

11

1

ˆ

=

−−−=∴

+

−=+

−

−−−

→

O

zO

O

K

zzzK

krKkk

HGIC

HxHEGxG43421

- 20 -

Fig. 16. Respuesta del sistema de control ante cambios de referencia y perturbación, utilizando un vector óptimo, E = [0.0079 0.0123]T y observadores de estados.

Vector óptimo 3

[ ]

×=

−−+

=

×=

×=

−

−−−−

0.0000 0.0000- 3.6099- 0.0000

0.0000 0 1.0000- 0

0.0000- 0.0000 0.0000 0.0000-

0.0000- 0.0000 0.0000 0.0000

101

1010

10100

010

10

11

10

10

H

AQA

ABBRQABBRAH

R

Q

TT

TTTT

Calculando los valores y vectores propios de H, tenemos:

6529.0

5.1

0002.0

104331.4

4

3

2

31

===

×=

λλλλ

=

0.6398 0.8641 1.0000- 0.9637

0.7686- 0.5034 0.0001 0.2670

0.0000 0.0000 0.0000- 0.0000-

0.0000- 0.0000- 0.0000 0.0000-

P

- 21 -

Considerando los valores propios cuya magnitud es menor que la unidad y sus correspondientes vectores, formamos la matriz T.

∴

−−−=

−−

−−

=

I

I

Y

X

T

6398.01

7686.00001.0

0000.00000.0

0000.00000.0

−−

=

−−

=

6398.01

7686.00001.0

0000.00000.0

0000.00000.0

I

I

Y

X

−−

==

=

=

4177.00002.0

5018.00000.0

6529.00

00002.0

0

0

4

3

H

H

EYW

E

I

λλ

11 −−−= IWXBRE T

[ ]0032.00021.0=E

( ) ( ) ( )

( ) ( )002901.0

1

1ˆ1lim

1

1

1

11

1

ˆ

=

−−−=∴

+

−=+

−

−−−

→

O

zO

O

K

zzzK

krKkk

HGIC

HxHEGxG43421

5.1. COMPARACIÓN DEL CONTROLADOR ÓPTIMO CON UNO NO ÓPTIM O Para comparar el desempeño del sistema de control con los vectores de ganancias óptimos, se presenta en la figura 18 la comparación de los tres sistemas calculados en la sección anterior con el desempeño de un sistema de control por realimentación de estados con el vector de ganancias E calculado utilizando la fórmula de Ackermann, y los polinomios normalizados de Bessel para un tss = 0.5 [s]. Fácilmente se puede comprobar que la respuesta del sistema de control con el vector óptimo 3, tiene un tiempo de establecimiento menor a 1 segundo y prácticamente la respuesta a las perturbaciones es inmediata, además de no ofrecer sobreimpulsos, tanto para cambios de referencia, como para perturbación. Las repuestas de los sistemas de control con los vectores óptimos 1, 2 priorizan la velocidad en la respuesta, sin embargo poseen un sobreimpulso mayor que la respuesta del sistema utilizando el vector calculado con la fórmula de Ackermann y el polinomio normalizado de Bessel para un tiempo de establecimiento, tss = 0.5 [s].

- 22 -

Se presenta, asimismo en la figura 19 las gráficas de comparación del estado X2 del sistema. Nuevamente se puede notar el mejor comportamiento del sistema de control con vector de ganancias óptimo 3.

Fig. 17. Respuesta del sistema de control ante cambios de referencia y perturbación, utilizando un vector óptimo, E = [0.0021 0.0032]T y observadores de estados.

Fig. 18. Comparación de la respuesta de los sistemas de control con diferentes vectores de realimentación, ante cambios de referencia y perturbación.

- 23 -

Fig. 19. Comparación de la respuesta del estado X2 de los sistemas de control con diferentes vectores de realimentación, ante cambios de referencia y perturbación.

Para la implementación práctica, se procedió de manera idéntica a la implementación del sistema de control por realimentación de estados, el único cambio radica en utilizar el nuevo vector de realimentación óptimo. En la figura 20 se ilustra el desempeño del sistema de control óptimo y en la figura 21, se presenta el desempeño del sistema de control utilizando el vector calculado por la fórmula de Ackermann y los polinomios de Bessel, tss = 0.5 [s].

Fig. 20. Desempeño del sistema de control óptimo en el control de velocidad de un motor c.c.

- 24 -

Fig. 21. Desempeño del sistema de control no óptimo en el control de velocidad de un motor c.c.

Claramente se puede notar del análisis de las figuras 20 y 21, el mejor desempeño del sistema de control óptimo, puesto que presenta una respuesta sin sobreimpulsos ante cambios de referencia y perturbación, además de estabilizarse más rápidamente que el sistema de control no óptimo. 6. CONTROL DE MODELO PREDICTIVO El control de modelo predictivo basa su principio de funcionamiento en el pronóstico del comportamiento de la planta para señales de control u[k] estimadas en un marco finito de optimización, a través del modelo de ésta, que como se podrá suponer tendrá un impacto directo en el desempeño del controlador. Los principios básicos del MPC son:

− A cada instante de muestreo evaluar la situación del sistema de control.

− Basado en esta evaluación, calcular la mejor estrategia.

− Aplicar solamente la primera muestra del cálculo de la estrategia, Receding Control.

− Re-evaluar la situación del sistema de control.

− Recalcular la mejor estrategia.

− Aplicar el criterio de Receding Control.

− Continuar con el proceso de manera progresiva. Es posible, por tanto suponer que son necesarios tres requisitos para establecer control MPC:

1) Estimación del comportamiento de la planta en un horizonte de predicción NP, a través de la estimación de la variable de control u[k] en un horizonte de control NC. Para éste paso, se requiere de un modelo preciso de la planta.

2) Evaluación del comportamiento real con respecto al deseado en la muestra Ki, de forma que sea posible recalcular la señal de control u[k] para alcanzar el comportamiento deseado. Se vuelve indispensable para éste objetivo poder medir con precisión la señal de respuesta del sistema a través de una etapa instrumental adecuada.

3) Establecer una estrategia de control realizable. Para éste diseño se utiliza una etapa integradora para contrarrestar el error de estado estable en la planta y el criterio de optimización en el horizonte de predicción se reduce a la minimización del índice de desempeño J = (Rs - Y)T(Rs - Y) + ∆UTR∆U.

- 25 -

Se presentan a continuación los cálculos del controlador MPC para la planta que identificamos en la sección 1, cuyo modelo se repite en la ecuación (12):

( )( ) 344798.02447.1

5831.312929.262

21

+−+= −

zz

zzz

zU

zY

( ) ( ) ( )( ) ( )

[ ]

( )( )

( )( ) ( )

( ) [ ] ( )( )

[ ]

( ) ( )[ ]

=

=

−+=∆

=

=

−=

∆=

∆

+

∆

=

++∆

=

=

−=

=+=+

−−−−

−

BCABCABCABCA

CBABCBCA

CBCAB

CB

Φ

CA

CA

CA

F

FxRΦRΦΦU

CBA

x0

BC

Bx

1AC

0Ax

CBA

xC

BxAx

BA

CPPPP

P

NNNNNN

isTT

mm

mm

mm

mm

Tmmm

mmm

mm

mmmm

k

ky

kky

kuky

k

ky

k

kky

kukk

121

22

1

0

0

00

000

100

29.26

58.31

29.26

112447.1

003448.0

012447.1

1

1

1

0158.31

29.26

03448.0

12447.1

1

OMMM

L

L

L

M

434214434421

×=

×=

×=

×=

1.5029

1.9135

2.3528

2.8159

101

101

0.1503 1.0995 1.0178

0.1914 1.3539 1.2576

0.2353 1.6008 1.4925

0.2816 1.8300 1.7133

101

4.2210 4.9496 5.5751 6.1024

4.9496 5.9074 6.7536 7.4767

5.5751 6.7536 7.8372 8.7879

6.1024 7.4767 8.7879 9.9816

101

3

6

4

5

sT

w

T

T

r

RΦ

FΦ

ΦΦ

(12)

- 26 -

La ley de control está dada por la siguiente ecuación:

( ) ( ) ( )11 +∆+=+ kukuku

En la figura 22 se ilustra el comportamiento del sistema de control utilizando la ley de control (13) y los cálculos bosquejados anteriormente. Resulta sencillo notar el excelente desempeño del controlador puesto que su tiempo de estabilización es menor a 1 segundo con un reducido sobreimpulso y su respuesta ante perturbaciones es menor a 1 segundo también. Claramente supera el desempeño de funcionamiento a los controladores implementados anteriormente: Realimentación de estados y control óptimo.

Fig. 22. Simulación del sistema de control de velocidad con MPC ante cambios de referencia y perturbación.

Utilizando LabVIEW como plataforma de desarrollo, los cálculos bosquejados en esta sección y la ley de control representada en la ecuación (13), se implementó el sistema de control MPC para regular la velocidad de la flecha del motor de corriente continua identificado en la sección 1. La curva de respuesta del sistema de control se ilustra en la figura 23. Resulta sencillo notar del análisis de la figura 23, que el comportamiento del sistema de control es muy parecido al presentado en la simulación de la figura 22. La respuesta ante cambios de referencia y perturbación es claramente superior a la de los otros sistemas de control implementados anteriormente, que se ilustran en las figuras 8, 20 y 21. El error de estable es “increíblemente” bajo reportándose mediciones de velocidad de 2000.4 rpm cuando la referencia es 2000 rpm, lo cual por supuesto no sucede en otros sistemas de control. El error de estado estable es de 2 x 10-4 %.

(13)

- 27 -

Fig. 23. Respuesta del sistema de control de velocidad con MPC ante cambios de referencia y perturbación.

7. SISTEMA DE DETECCIÓN DE FALLAS EN UN MOTOR C.C.

Los sistemas dinámicos utilizados en procesos productivos, por ejemplo, están sometidos a diversas fallas debidas a muchos factores tales como: mal uso, antigüedad, falta de mantenimiento, etc. En el caso particular de un motor, que se utiliza como elemento motriz en la industria, automóviles, manufactura, etc., puede presentar diferentes fallas debidas a causas variadas, que van desde alimentación inadecuada, desgaste de rodamientos, sobrecargas en el eje, entre otros. Resulta, por tanto muy importante desarrollar un sistema de detección de fallas confiable, y es justamente allí en donde enfocamos el desarrollo del presente trabajo, a través del uso de observadores de estado para fallas bien conocidas en el sistema. Todos los modelos aquí presentados se obtuvieron utilizando el método de identificación de mínimos cuadrados, que se explica en detalle en la sección 1. Repitiendo el procedimiento para los casos en los que el sistema presenta fallas, los modelos fueron los siguientes:

− Carga tipo escalón en el eje del motor. Mientras el motor se encontraba trabajando a una velocidad estable, se aplicó primeramente una carga pequeña al eje del motor y se utilizó el proceso de identificación. Se realizaron dos modelaciones, considerando en una segunda oportunidad una carga fuerte al eje del motor, generándose los siguientes modelos:

21

21

2

21

3

21

21

2

21

1

654276.057142.11

0224.1306579.0

654276.057142.1

0224.1306579.0)(

448921.033245.11

9834.376978.9

448921.033245.1

9834.376978.9)(

−−

−−−

−−

−−−

+−+=

+−+=

+−+=

+−+=

zz

zz

zz

zzzzHGf

zz

zz

zz

zzzzHGf

− Perturbación en el sensor de velocidad. En el encoder del motor se procedió a obstruir algunas ranuras del disco, cambiando ligeramente la linealidad de la medición y también la velocidad calculada. Se procedió, asimismo a aplicar la metodología de identificación, obteniendo el modelo siguiente:

21

21

2

21

2 561823.03951.11

618.142715.18

561823.03951.1

618.142715.18)( −−

−−−

+−+=

+−+=

zz

zz

zz

zzzzHGf

- 28 -

En virtud de que el objetivo principal es identificar las fallas que modelamos en el motor, el diseño en base al cual se realizan los cálculos y la implementación es el que se ilustra en la figura 22. Analizando la figura 22, podemos notar que el “secreto” del sistema se encuentra en el banco de observadores que estiman los valores de los estados para cada modelo obtenido en el punto anterior, y en el caso de presentarse una de las fallas modeladas, automáticamente el observador correspondiente convergerá a la variable medida, que en este caso es la velocidad de la flecha del motor, permitiéndonos emitir un diagnóstico y/o alarma de lo que ocurre en el sistema.

Figura 24. Banco de observadores para el sistema de detección de fallas.

− HGf1

[ ]

[ ]

==

∆=

=

+

−−

=

++

−

133245.1

01

)(

)(

)(01)(

)(9834.37

6978.9

)(

)(

0448921.0

133245.1

)1(

)1(

1

2

1

2

1

2

1

TCGCWo

JnWoGcKe

kx

kxky

kukx

kx

kx

kx

M

- 29 -

Considerando Deadbeat en la ubicación de polos:

−=

= −

448921.0

33245.12

12

Ke

JWoGKe

− HGf3

[ ]

[ ]

==

∆=

=

+

−=

++

−

157142.1

01

)(

)(

)(01)(

)(0224.13

06579.0

)(

)(

065427.0

157142.1

)1(

)1(

1

2

1

2

1

2

1

TCGCWo

JnWoGcKe

kx

kxky

kukx

kx

kx

kx

M

Considerando Deadbeat en la ubicación de polos:

−=

= −

65427.0

57142.12

12

Ke

JWoGKe

− HGf2

[ ]

[ ]

==

∆=

=

+

−=

++

−

139531.1

01

)(

)(

)(01)(

)(618.14

2715.18

)(

)(

0561823.0

139531.1

)1(

)1(

1

2

1

2

1

2

1

TCGCWo

JnWoGcKe

kx

kxky

kukx

kx

kx

kx

M

Considerando Deadbeat en la ubicación de polos:

−=

= −

561323.0

39531.12

12

Ke

JWoGKe

En la figura 24, se presenta la simulación del observador de estados de la planta. Claramente se puede concluir que el observador calcula de manera correcta los estados de la planta, puesto que converge al resultado ideal.

- 30 -

Fig. 25. Desempeño del observador de estados de la planta.

El desempeño de los observadores de los modelos detallados previamente, se presentan en las figuras 25, 26 y 26. Cabe destacar que en éstas gráficas la planta en principio no se encuentra sometida a ninguna perturbación y en el tiempo 2.5 [s] se presenta la falla, permitiendo que el observador respectivo converja.

Fig. 26. Desempeño del observador de estados de la planta + falla 1 (carga leve en el eje del motor).

- 31 -

Fig. 27. Desempeño del observador de estados de la planta + falla 3 (carga moderada en el eje del motor).

Fig. 28. Desempeño del observador de estados de la planta + falla 2 (perturbación en el sensor).

El comportamiento de todos los observadores para los diferentes casos: sin falla, falla 1, falla 2 y falla 3, se presenta en las figuras 29, 30 y 31. Resulta sencillo notar allí que para el caso concreto de la falla presente, el observador correspondiente converge y los demás no, siendo sencillo a la larga a través de comparaciones entre los residuos de los observadores y la velocidad del motor determinar que falla se encuentra activa en el sistema.

- 32 -

Fig. 29. Desempeño de los observadores de estados de la planta + falla 1 (carga leve en el eje del motor).

Fig. 30. Desempeño de los observadores de estados de la planta + falla 2 (perturbación en el sensor).

- 33 -

Fig. 31. Desempeño de los observadores de estados de la planta + falla 3 (carga moderada en el eje del motor).

Utilizando LabVIEW como plataforma de desarrollo, se implementaron los observadores calculados en la sección 3, obteniéndose los resultados presentados en las figuras 32, 33, 34 y 35. Claramente se puede notar del análisis de éstas gráficas que los resultados de las simulaciones concuerdan con los de la práctica, validando de ésta manera el diseño y los cálculos del banco de observadores planteado en la figura 3.

Figura 32. Desempeño de los observadores de estados en ausencia de fallas.

- 34 -

Figura 33. Desempeño de los observadores de estados ante la falla 1.

Figura 34. Desempeño de los observadores de estados ante la falla 2.

- 35 -

Figura 35. Desempeño de los observadores de estados ante la falla 3.

8. CONCLUSIONES La elaboración de este proyecto de control trajo algunas conclusiones que son de gran utilidad para el diseño y aplicación de los sistemas de control en algún proceso dado. A lo largo de este documento se presentaron los diseños de varios controladores para una misma planta, esto es de gran beneficio, ya que se puede hacer una comparación del desempeño de los diferentes controladores bajo las mismas condiciones de operación, su complejidad de diseño y por tanto de una forma indirecta, su margen costo-beneficio para una aplicación en especial. La primer parte del proyecto se refirió a la identificación de la planta, esto es, la obtención del modelo matemático de 2do orden que describiera la dinámica del motor de corriente directa que se estuvo usando. Esta etapa del proyecto es de vital importancia, ya que si un modelo matemático que se obtenga de la planta no describe el funcionamiento real de la misma, ningún algoritmo de control diseñado basado en ese modelo matemático será de utilidad para tener el control de la planta. El método utilizado para hacer la identificación del modelo matemático fue el de “Mínimos Cuadrados por Lotes” con datos entrada salida proporcionados por la aplicación de una señal PRBS alrededor de punto de operación de 2000 rpm, lo que proporcionó un modelo de manera sencilla. El primer controlador desarrollado fue el de realimentación de estados. Considerando que los estados no siempre están disponibles para su medición, se propuso el diseño de observadores de estados (Obtenidos mediante la fórmula de Ackermann). De acuerdo a la literatura citada, el tiempo de estabilización para el diseño de los observadores de estado debe de ser entre 4 y 6 veces más rápida que la respuesta del sistema de control, esto para evitar hacer que la dinámica de respuesta del observador tenga dominio sobre la respuesta general del sistema.

- 36 -

Al sistema de retroalimentación de estados se le agrego un sistema de seguimiento para que el sistema pudiera seguir las referencias arbitrarias sin error de estado estable. Sin embargo, esta dinámica adicional que se agregaba para hacer el seguimiento producía una respuesta no tan favorable cuando se presentaban perturbaciones en el sistema. Después de desarrollar el controlador mediante la técnica de retroalimentación de estados y el uso de observadores, se propuso un sistema que identificara fallas presentes en la planta. Estas fallas producían una cierta dinámica en la planta, con lo que el modelo del sistema original se modificaba. Al conocer de manera matemática la dinámica de las fallas presentes en el sistema, se creó un banco de observadores con las dinámicas de cada una de las fallas que podrían estar presentes. Al inducir cada una de las fallas mientras el sistema se encontraba trabajando, se pudo comprobar mediante la simulación y mediante la implementación física como era que el banco de observadores ayudó a identificar la dinámica que estaba presente, ya fuera la del sistema trabajando en forma normal y la de la presencia de alguna de las fallas inducidas. Cabe destacar que al igual que con la primer identificación del sistema, la correcta modelación de la dinámica del sistema cuando una falla estaba presente fue necesaria, debido a que si la el sistema presentaba una falla por alguna razón en particular y esa falla no está bien descrita en la dinámica del modelo adquirido, el sistema de detección de fallas nunca se daría cuenta que esa falla estaba presente por el error presente en el “observador” de esa falla en particular. Adicionalmente a los sistemas anteriores, se propuso el diseño de un controlador óptimo cuadrático, mediante la técnica de minimización al utilizar los multiplicadores de Lagrange. Un controlador óptimo está enfocado a minimizar o maximizar un cierto índice de desempeño en particular. Se debe de considerar que al optimizar un cierto índice de desempeño, no siempre optimiza otro índice de desempeño, esto es por ejemplo, minimizar el tiempo de respuesta de un sistema regularmente no siempre será el óptimo para la vida útil de actuador que hace que el tiempo de respuesta del sistema sea rápido. Una conclusión importante en este tipo de controlador, es que al trabajar con una ley de control optimo basada en un índice de desempeño que se basa en un índice de desempeño cuadrático en la cual se obtiene como resultado una función lineal del vector de estados x(k). Si se hace una comparación entre el desempeño del controlador optimo con el no optimo, se puede llegar a conclusiones interesantes. En el controlador óptimo se puede ver un desempeño muy suave en comparación con el desempeño provocado por los polinomios de Bessel (no óptimo). En la respuesta del controlador óptimo prácticamente no existe sobretiro y la respuesta ante perturbaciones es prácticamente instantánea, en comparación con el control no óptimo, aunque también se vio que si el índice que se desea minimizar es el tiempo de respuesta, el vector de retroalimentación óptimo produjo una respuesta con sobretiro mayor a la respuesta del sistema no óptimo. El control óptimo pretende minimizar la diferencia del comportamiento del sistema real, para lo cual la propuesta de matrices Q y R, en el caso de control optimo cuadrático, se vuelve el factor clave de diseño. En el caso de que se requiera penalizar amplitudes de un estado en específico, se debe recurrir a incrementar el índice correspondiente de la matriz Hermítica definida positiva Q, considerablemente. En el desarrollo de éste trabajo, se utilizaron incluso cantidades de hasta 1 × 10

��. Los requisitos básicos para la implementación del MPC son: modelo, medición y ley de control realizable. En el caso de nuestro proyecto el modelo está basado en una estructura ARX, la medición se realiza con un encoder incremental instalado en el eje del motor y la ley de control está basada en el criterio de optimización

en el horizonte de predicción minimizando el índice de desempeño J = (Rs - Y)T(Rs - Y) + ∆UTR∆U.

- 37 -

9. BIBLIOGRAFÍA

[1] Ioan D. Landau and Gianluca Zito. Digital Control Systems. Lab. d’Automatique de Grenoble. Springer, London – 2006.

[2] Antonio Favela C., Método de Identificación por Mínimos Cuadrados, ITESM – Cátedra de

Identificación de Procesos.

[3] Ogata Katsuhiko, Sistemas de control en tiempo discreto, segunda edición, Prentice – Hall, 1996. [4] Santina, Stubberud, Hostetter , Digital Control Systems Design, 2nd. Edition, Saunders College

Publishing, 1994.

[5] Morillo, Ríos-Bolíva, Control of a DC Brush Motor via Dynamic Output Feedback, Zulia University, Venezuela.

[6] Mosca Edoardo, Optimal, Predictive and Adaptive Control, Prentice-Hall, 1995

[7] Wand Liuping, Model Predictive Control System Design and Implementation Using Matlab, Springer,

Australia – 2009.