ANÁLISIS DE CONTROLABILIDAD PARA UN MODELO DINÁMICO DE UNA …

1

ANÁLISIS DE CONTROLABILIDAD PARA UN MODELO DINÁMICO DE UNA

MOTOCICLETA

Trabajo de grado Nº 0966

CARLOS ANDRÉS QUIROGA RENGIFO

PROYECTO DE GRADO PRESENTADO PARA OPTAR AL TÍTULO DE INGENIERO

ELECTRÓNICO

DIRECTOR

INGENIERO KAMILO ANDRES MELO BECERRA M.Sc

PONTIFICIA UNIVERSIDAD JAVERIANA

FACULTAD DE INGENIERÍA

CARRERA DE INGENIERIA ELECTRÓNICA

BOGOTÁ D.C.

MAYO DE 2010

PONTIFICIA UNIVERSIDAD JAVERIANA

2

FACULTAD DE INGENIERÍA

CARRERA DE INGENIERIA ELECTRÓNICA

Rector de la Universidad:

Padre Joaquín Sánchez García S.J.

Decano Académico de la Facultad de Ingeniería:

Ingeniero Francisco Javier Rebolledo Muñoz

Decano Medio Universitario Facultad de Ingeniería:

Padre Sergio Bernal S.J.

Director de la Carrera de Ingeniería Electrónica:

Ingeniero Juan Manuel Cruz Bohórquez

Director de Departamento de Ingeniería Electrónica:

Ingeniero Jorge Luis Sánchez Téllez

3

Artículo 23 de la resolución No. 13 del 6 de julio de 1964:

“La universidad no se hace responsable de lo conceptos emitidos por sus alumnos en sus proyectos de grado.

Solo velará porque no se publique nada contrario al dogma y la moral católica y porque no contenga ataques o polémicas puramente personales. Antes bien, que se vean en ellos el anhelo de buscar la verdad y la justicia.”

4

TABLA DE CONTENIDO 1. INTRODUCCIÓN ................................................................................................................................. 5

2. MARCO TEÓRICO ............................................................................................................................... 6

2.1. Breve historia de la motocicleta y de sus primeros modelos matemáticos .................................... 6

2.2. Ecuación de Euler-Lagrange y el Lagrangiano .............................................................................. 8

2.3. Modelo de punto de masa y sus ecuaciones dinámicas .................................................................. 9

2.4. Análisis de la controlabilidad de un sistema MIMO ....................................................................13

3. DESARROLLOS .................................................................................................................................15

3.1. Construcción del modelo en Simulink™ .....................................................................................15

3.1.1. Planteamiento de las ecuaciones para la construcción del modelo ......................................15

3.1.2. Simplificación del modelo ...................................................................................................18

3.1.3. Adaptación del modelo .........................................................................................................19

3.1.4. Validación del modelo simplificado .....................................................................................20

3.2. Linealización ................................................................................................................................28

3.3. Controlabilidad .............................................................................................................................30

3.3.1. Escalización ..........................................................................................................................31

3.3.2. Realización Mínima .............................................................................................................32

3.3.3. Controlabilidad funcional .....................................................................................................33

3.3.4. Obtención de los polos del sistema ......................................................................................35

3.3.5. Obtención de los ceros del sistema .......................................................................................37

3.3.6. Matriz RGA y su respuesta en frecuencia ............................................................................38

3.3.7. SVD y respuesta en frecuencia de los valores singulares .....................................................40

4. CONTROLADOR PROPUESTO ........................................................................................................42

5. ANÁLISIS DE RESULTADOS Y CONCLUSIONES .......................................................................45

6. BIBLIOGRAFÍA ..................................................................................................................................48

7. ANEXOS ..............................................................................................................................................49

5

1. INTRODUCCIÓN

El presente es un proyecto que se enmarca en el área de control y se enfoca en el tema específico de las motocicletas. Está fuertemente influenciado por [1], el libro de control multivariable escrito por Sigurd Skogestad e Ian Postlethwaite. Pretende continuar y complementar los múltiples estudios sobre motocicletas que se vienen adelantando en otros países [2] [3] y a la vez motivar en Colombia el desarrollo de aplicaciones afines (A este respecto cabe anotar que en la Pontificia Universidad Javeriana se han realizado recientemente tres trabajos de grado aplicados a las motocicletas). Está relacionado con los sistemas dinámicos multivariados y tiene grandes aplicaciones a futuro en las áreas de seguridad para motociclistas, vehículos no piloteados de exploración extraterrestre y sistemas de locomoción terrestre, donde el desempeño de un sistema de cuatro ruedas no es tan eficiente como el de uno de dos (por efectos de fricción, peso, gasto de combustible, etc.), entre otros.

Específicamente se va a realizar un análisis de controlabilidad para una motocicleta que se va a representar a partir de su modelo dinámico como un sistema MIMO de dos entradas y tres salidas. La idea es determinar si es posible mantener el equilibrio de este vehículo controlándolo mediante el torque del motor y el ángulo de su timón.

La metodología será recoger un modelo de la literatura y adaptarlo a las necesidades del proyecto. Se corregirán posibles imprecisiones limitando el modelo de acuerdo a las características físicas del sistema que se está representando y se modificará su dinámica, de ser necesario, para representar adecuadamente la respuesta del sistema real. Este modelo se construirá en Simulink™, un entorno de programación gráfica que hace parte del programa MATLAB®. También se tendrá que linealizar para luego analizar su controlabilidad, al término de lo cual se espera obtener información valiosa para proponer futuros desarrollos como la síntesis de un controlador o el análisis de robustez del sistema.

6

2. MARCO TEÓRICO

2.1. Breve historia de la motocicleta y de sus primeros modelos matemáticos

Los primeros vehículos impulsados por humanos fueron de cuatro ruedas. Se accionaban directamente por mecanismos de palancas y engranajes. La dificultad e incomodidad a la hora de usarlos le abrieron paso a una nueva generación de sistemas de transporte impulsados por humanos. Estos nuevos sistemas contaban con sólo dos ruedas.

El primer diseño lo desarrolló el inventor alemán Karl von Drais en 1816 y se conocía como “máquina de correr” o Draisine (ver Figura 1). Su construcción era muy similar a la de la bicicleta de la actualidad. El marco y las ruedas eran de madera, tenía silla y descansabrazos, la rueda delantera podía girar para cambiar la dirección de la máquina y no tenía pedales por lo que el piloto se debía impulsar como si estuviera corriendo. Dentro de sus fallas más notorias estaba la falta de frenos y de un mejor sistema de propulsión.

Figura 1 Tomado de [9]

Al diseño del Barón Karl von Drais se le hicieron bastantes mejoras basadas tanto en la experiencia con la máquina como en las nuevas tecnologías y materiales que se fueron desarrollando. Dentro de los cambios más significativos se encuentran la implementación de pedales con cadenas y engranajes como sistema de propulsión, el diseño e implementación de ruedas neumáticas y marcos de metal más livianos.


En 1867 el inventor estadounidense Sylvester Howard Roper le añadió un motor a vapor al denominado “bone shaker” (uno de los diseños a los que se llegó a partir del Draisine) creando lo que se podría considerar como la primera motocicleta (ver Figura 2). Esta era una máquina impulsada por un motor a vapor de dos cilindros, alimentado con carbón, que accionaba una manivela, conectada a la rueda trasera, a través de dos varillas. Sin embargo, Gottlieb Daimler es considerado por muchos como el inventor de la primera motocicleta porque su diseño fue el primero en emplear un motor de combustión interna para su

7

propulsión (ver figura 3). Este desarrollo se dio en 1885. Sus llantas estaban construidas en madera (al igual que el marco) recubiertas por una banda de hierro, también tenía llantas laterales en madera.


Con el paso del tiempo, las nuevas tecnologías y las necesidades de los nuevos usuarios se desarrollaron diseños más especializados para, por ejemplo y entre otras cosas, alcanzar una mayor velocidad a un menor costo. Esto se fue logrando con motores más eficientes, diseños con mayor estabilidad, etc. Aún en la actualidad se sigue intentando mejorar este medio de transporte. Desde el punto de vista del control, y teniendo en cuenta que cada vez son menos evidentes las posibles mejoras, se hace necesario contar con modelos que representen el sistema adecuadamente, una teoría de control sólida y un análisis dinámico riguroso para lograr nuevos adelantos.

En 1869 se publicaron cinco artículos sobre las dinámicas de una bicicleta. Estos artículos constituyen el análisis más antiguo que se conoce sobre el tema y están basados en el modelo heurístico de un péndulo invertido. A través de este modelo el autor pretendió estudiar el balance, la dirección y la propulsión de la bicicleta, sin embargo hoy en día estos artículos no tienen mucho valor desde el punto de vista técnico. En cuanto a aportes técnicos se refiere, la primera contribución importante estuvo a cargo del matemático Francis John Welsh Whipple en 1899. Este trabajo incluye por primera vez un conjunto de ecuaciones diferenciales no lineales que describen el movimiento de la bicicleta y su piloto. Whipple modeló la bicicleta como dos marcos unidos por una bisagra. Las llantas trasera y delantera hacen parte de los marcos trasero y delantero respectivamente y cada uno puede rotar libremente con respecto del otro. El piloto es una masa inerte unida al marco trasero. Se asume que las llantas no se resbalan y que son tan delgadas que tocan el piso en un solo punto (cada una). Estas suposiciones se siguen usando en modelos más completos que están vigentes en la actualidad.

Después del modelo descrito en el párrafo anterior se propusieron otros. Algunos de estos incluían elementos tomados del trabajo de Whipple para modelar las dinámicas del vehículo mientras otros se diseñaron para modelar el comportamiento de partes específicas de este, como las llantas o la suspensión. Dadas las similitudes que existen entre las bicicletas y las motocicletas, el modelo de Whipple se puede adaptar para analizar las dinámicas de estas últimas. De hecho en este trabajo se va a usar el modelo del punto de masa, que se puede ver como un caso especial del modelo de Whipple. Este modelo es ampliamente empleado para modelar motocicletas y se ha comprobado que, a pesar de todas las suposiciones que se hacen, los resultados que se obtienen al usarlo son muy cercanos a la realidad. Más adelante se presentará formalmente el modelo que se utilizará a lo largo de este trabajo.

8

2.2. Ecuación de Euler-Lagrange y el Lagrangiano

En esta sección, y debido a la importancia de este concepto en la deducción de las ecuaciones dinámicas que describen el sistema propuesto, se pretende ilustrar el uso del Lagrangiano y cómo se puede obtener a partir de la segunda ley de Newton. También se muestra cómo se obtiene la ecuación de Euler-Lagrange. Para esto se toma un ejemplo sencillo que aparece en [4]. Allí mismo se encuentran una explicación más profunda y ejemplos más complejos en caso tal que el lector quiera complementar la información que se presenta a continuación.

Figura 4 – Diagrama de fuerzas de un sistema con un grado de libertad compuesto por una partícula de masa m que se mueve sobre el eje y.

Considérese un sistema con un grado de libertad compuesto por una partícula de masa m que se mueve sobre el eje y. Sobre la partícula actúan una fuerza f en el sentido positivo del eje mencionado y la fuerza de gravedad mg en el sentido opuesto (ver Figura 4). De acuerdo con la segunda ley de Newton, se obtiene la siguiente ecuación:

�� = � − ��

Ecuación 1

La parte izquierda de la ecuación se puede escribir como sigue:

�� = �� = �� 12 �� = �� Ecuación 2

En la Ecuación 2 se ha usado k para representar la energía cinética. Además, teniendo en cuenta que en la gran mayoría de sistemas la energía cinética depende de más de una variable, se usó la notación de derivada parcial.

También la fuerza de gravedad se puede escribir de la siguiente forma:

�� = � �� = � �

Ecuación 3

Donde P representa la energía potencial de la partícula.

9

Teniendo los términos de la ecuación escritos como las energías cinética y potencial del sistema es posible escribir el Lagrangiano del mismo, pues éste se define como la diferencia entre la energía cinética y la energía potencial (Ecuación 4).

� = � − � = 12 �� − ��

Ecuación 4

Para obtener las igualdades que se presentan en la Ecuación 5 se deriva la Ecuación 4 con respecto a la primera derivada de y. Luego se repite la operación, pero derivando con respecto a y. Al observar estas derivadas y compararlas con la Ecuación 2 y con la Ecuación 3 se llega a lo siguiente:

� �� = � �� y � � = � �

Ecuación 5

Reemplazando los términos que conforman la Ecuación 1 por los términos equivalentes que se muestran en la Ecuación 5 y reorganizándola, se obtiene lo que se conoce como la ecuación de Euler-Lagrange del sistema. A continuación se presenta la ecuación de Euler-Lagrange que se obtuvo para este sistema:

�� − � � = �

Ecuación 6

Aunque esta ecuación es equivalente a la ecuación que describe la segunda ley de Newton, cuando se trabaja con sistemas que dependen de más de una variables tiene ventajas operativas sobre las ecuaciones vectoriales de Newton. Una de las ventajas es que al aplicar esta ecuación para cada coordenada generalizada se obtiene un sistema de ecuaciones diferenciales cuya solución expresa todas las fuerzas generalizadas que interactúan. Adicionalmente la conservación del momento generalizado se puede verificar observando las dependencias de la función de energía cinética, pues si una función tiene dependencia de las velocidades generalizadas, pero no de la posición existirá conservación del momento.

En este caso se llegó a la ecuación de Euler-Lagrange, pero no siempre es así. Un ejemplo es la motocicleta que se estudiará en este proyecto. Por ser un sistema no holonómico tiene ciertas restricciones. Al tener en cuenta las restricciones no holonómicas de un sistema, su Lagrangiano se conoce como “Lagrangiano restringido” y las ecuaciones que se obtienen a partir de éste son equivalentes a las ecuaciones d’Alambert para sistemas restringidos y no a las de Euler-Lagrange, aunque son muy similares. Para profundizar en este punto se recomienda consultar [5].

2.3. Modelo de punto de masa y sus ecuaciones dinámicas

A continuación se describirá el modelo de punto de masa. Se empezará por enumerar las restricciones que se hicieron sobre el sistema para construir el modelo.

1. El sistema es no holonómico. 2. No existe ningún tipo de deslizamiento (ni derrape lateral ni deslizamiento rotacional) ni en la

llanta trasera ni en la delantera. 3. Las llantas tocan el suelo en un solo punto cada una, su grosor es despreciable.

10

4. La masa de la motocicleta está concentrada en un solo punto; no se consideran los momentos de inercia de las ruedas ni de ninguna otra parte de la motocicleta.

5. La motocicleta se desplaza sobre una superficie plana

Para empezar se considera un marco de referencia inercial con los ejes X y Y separados 90 grados entre sí generando el plano del suelo. Perpendicular a este plano y en dirección opuesta a la fuerza de gravedad se tiene el eje Z. La intersección del plano de simetría de la motocicleta con el suelo forma lo que se conoce como línea de contacto. Esta línea de contacto gira alrededor del eje Z creando un ángulo Ψ con respecto al eje X. Se asume que el ángulo Ψ es cero cuando la línea de contacto se encuentra paralela al eje X. El ángulo que se genera entre el eje Z y el plano de simetría de la motocicleta se denomina ángulo de inclinación y se representará con la variable Θ. A continuación se definen estas y otras variables antes de continuar con la presente descripción:

m Masa de la motocicleta

L Distancia entre los dos puntos de contacto de las ruedas con el suelo

h Altura del centro de masa de la motocicleta

b Distancia horizontal entre el punto de contacto de la rueda trasera y la proyección del centro de masa sobre el suelo

R Radio de giro de la motocicleta

η Ángulo de el tenedor delantero con respecto a la horizontal

∆ Distancia entre el punto de contacto de la rueda delantera y la proyección del tenedor delantero sobre el suelo

Θ Ángulo de inclinación de la rueda trasera de la motocicleta con respecto a la vertical

Ψ Ángulo de la rueda trasera de la motocicleta con respecto a la horizontal

Φ Ángulo de giro del timón

β Ángulo de giro de la rueda delantera (ángulo efectivo del timón)

En las Figuras 5 a) y b) se puede apreciar gráficamente lo que cada una representa.

Figura 5 a) Tomado de [3] Figura 5 b) Tomado de [3]

11

Continuando con la descripción del modelo, se tiene en cuenta también el ángulo β que se genera entre la intersección del plano de la llanta delantera y la línea de contacto. Con base en la parametrización de este ángulo se define una variable σ de la siguiente manera:

�: = 1� = tan ��

Ecuación 7

Se puede igualar la parte del extremo izquierdo de la Ecuación 7 y la parte del extremo derecho de esta

misma ecuación. Si se deriva esta igualdad con respecto al tiempo y se despeja �� se obtiene lo siguiente:

�� = �1 + tan� � �� = �1 + �� Ecuación 8

De aquí en adelante σ será denominada “variable de giro”. En este punto se tiene que las coordenadas generales de la motocicleta son x, y, Ψ, Θ y σ, siendo (x, y) las coordenadas cartesianas del punto de contacto de la llanta trasera con el suelo.

A partir de la geometría del mecanismo de giro de la llanta delantera se puede obtener la siguiente relación entre los ángulos Φ, β, Θ y η:

tan � cos # = tan ∅ sin &

Ecuación 9

Por otro lado, de la Figura 5 b) se pueden calcular las componentes en x y y de la velocidad lineal vr de la rueda trasera y de la velocidad lateral o perpendicular vp, de la misma, obteniendo la ecuación que se presenta a continuación:

'(�� ) = *cos + −sin +sin + cos + , *-.-/,

Ecuación 10

Por otro lado, teniendo en cuenta la restricción número 2 del modelo se deduce que la velocidad lateral o perpendicular -/ es igual a cero. La otra igualdad que aparece a continuación se deriva de las propiedades

no holonómicas del sistema (la motocicleta no puede cambiar su orientación sin desplazarse, por lo que la velocidad con la que cambie el ángulo + del sistema dependerá de la velocidad lineal -. de la motocicleta) y se obtiene tras inspeccionar la Figura 5 b).

0+� = -. tan �� = �-.-/ = 0 2 Ecuación 11

También se debe tener en cuenta el mecanismo de estabilización de la motocicleta a bajas velocidades mediante el “timoneo”, es decir el efecto del cambio en el ángulo Φ sobre el ángulo Θ. Para esto se puede calcular el cambio de altura del centro de masa de la motocicleta debido al movimiento del timón a partir

12

de una rotación del marco un ángulo δ en el plano del suelo XOY. En [12] se obtiene un aproximado de δ de la siguiente manera:

3 = ∆ sin &� �

Ecuación 12

De acuerdo con esto último se puede calcular la reducción en la altura del centro de masa como en la Ecuación 13.

ℎ6∆ = 37 sin # = 7∆ sin &� � sin # ≈ 7∆� sin & sin #

Ecuación 13

Con base en todo lo anterior se puede deducir el Lagrangiano resultante del presente sistema como:

�9 = −��ℎ cos # − 7∆� sin & sin #� + 12 :;∅� �+ 12 � *-. + ℎ�-. sin #�� + ℎ�#� � sin� # + <7�-. − ℎ#� cos #=�,

Ecuación 14

En esta ecuación Js es el momento de inercia del mecanismo de giro.

A partir del Lagrangiano obtenido se pueden deducir las ecuaciones que describen la dinámica de la motocicleta de la siguiente manera:

�� = >?

Ecuación 15

M ' #�-�.) = K + B *>?C. ,

Ecuación 16

Como se puede ver ωσ es la entrada virtual del sistema que representa la velocidad angular del timón y Fr

es la fuerza de tracción longitudinal de la llanta trasera. Las demás componentes de la ecuación son:

M = ' ℎ� −7ℎ� cos #−7ℎ� cos # 7�� + 1 + ℎ� sin #��)

Ecuación 17

B = D 7ℎ-. cos # 0−E7�� + ℎ sin #1 + ℎ� sin #�F-. 1�G

Ecuación 18

K = H�ℎ sin # + 7∆� sin & cos #� + 1 + ℎ� sin #�ℎ�-.� cos #−2ℎ�-.#�1 + ℎ� sin #� cos # − 7ℎ�#� � sin # I

Ecuación 19

13

2.4. Análisis de la controlabilidad de un sistema MIMO

En esta sección se presenta un procedimiento tomado de [1] para el análisis de la controlabilidad de un sistema MIMO. Este proceso presupone que se han escogido las entradas y salidas de la planta y se quiere analizar su modelo G para determinar qué desempeño en términos de control se puede esperar de esta.

1. Escalizar todas las variables (entradas u, salidas y, disturbios d, referencias r) para obtener el modelo escalizado, y = G(s)u + Gd(s)d, r = Rř, donde G(s), Gd(s) y r son el modelo de la planta, el de los disturbios y la referencia respectivamente.

2. Obtener la realización mínima del sistema. 3. Revisar la controlabilidad funcional. Para poder controlar las salidas independientemente, lo

primero que se necesita es tener tantas entradas u como salidas y. Lo segundo es que el rango de G(s) sea igual al número de salidas l. Es decir, que el valor singular mínimo de G(jω), σ(G)

= σl(G), debe ser diferente de cero (excepto en los ceros de los ejes jω). Si la planta no es controlable funcionalmente, encuentre la dirección de salida en que la planta no tenga ganancia para obtener información sobre la causa del problema.

4. Localice los polos. Se debe obtener la ubicación de los polos RHP (inestables) y sus direcciones asociadas. Los polos “rápidos” alejados del origen son perjudiciales.

5. Localice los ceros. Se debe obtener la ubicación de los ceros RHP y sus direcciones asociadas. También es preciso buscar los ceros fijados en ciertas salidas. Los ceros “pequeños” (cercanos al origen) son perjudiciales.

6. Obtener la respuesta en frecuencia G(jω) y construir la matriz RGA (Λ = G X (G-1)T). Las plantas con grandes elementos RGA a frecuencias de transición son difíciles de controlar y se deben evitar. (El operador X representa el producto punto entre matrices.)

7. De este punto en adelante la escalización es crítica. Se deben computar los valores singulares de G(jω) y se deben graficar como una función de la frecuencia. Los vectores singulares asociados de entrada y de salida también deben considerarse.

8. El valor singular mínimo, σ(G(jω)), es una medida de controlabilidad particularmente útil. Por lo general debe ser tan grande como sea posible a frecuencias en las que se busca ejercer control. Si σ(G(jω)) < 1 no será posible generar cambios independientes de magnitud unitaria a la salida a partir de entradas de magnitud unitaria a una frecuencia ω determinada.

9. En cuanto a los disturbios y la necesidad de control, considere la matriz Gd. A frecuencias donde uno o más elementos son mayores que 1, se necesita control. Se obtiene una mayor información examinando un disturbio a la vez (las columnas gd de la matriz Gd). Para cada disturbio se requiere que S sea menor que 1/||gd||

2 en la dirección yd del disturbio, es decir que ||Syd||

2 ≤ 1/||gd||2. En consecuencia se requiere que al menos σ(S) ≤ 1/||gd||

2 y puede que también se requiera que σ(S) ≤ 1/||gd||

2. 10. Disturbios y saturación a la entrada:

Primer paso: Considere las magnitudes de entrada necesarias para un control perfecto computando los elementos de la matriz G-1

Gd. Si todos los elementos son menores que 1 a todas las frecuencias, se espera que la saturación a la entrada no sea un problema. Si algunos elementos de G-1

Gd son mayores que 1 quiere decir que el control perfecto (e = 1) no se podrá lograr a esas frecuencias. Sin embargo no se descarta la posibilidad de lograr el “control aceptable” (||e||2 < 1). Lo anterior se podrá poner a prueba en el segundo paso. Segundo paso: Considere los elementos de UH

Gd y asegúrese de que los elementos de la i-ésima fila sean menores que σi(G) + 1 a todas las frecuencias.

14

11. Ver los disturbios, los polos RHP, los ceros RHP y sus ubicaciones y direcciones asociadas para determinar si los requerimientos son compatibles.

12. Si el número de condición γ(G) es pequeño, se espera no tener problemas con la incertidumbre. Si los elementos RGA son grandes, es de esperarse una alta sensibilidad del sistema a las incertidumbres.

Es importante aclarar que no se incluirá análisis de disturbios en este trabajo, por lo que no se seguirán los puntos del 9 al 12 incluidos en el análisis de controlabilidad que se acaba de describir.

15

3. DESARROLLOS

A partir de este punto se presentará todo el proceso que se adelantó para cumplir con los objetivos propuestos.

3.1. Construcción del modelo en Simulink™

Simulink™ es una herramienta que permite construir modelos que representan sistemas variantes en el tiempo para su análisis y simulación. Esta aplicación hace parte del programa MATLAB® y tiene un enfoque gráfico basado en diagramas de bloques donde se construyen modelos complejos a partir de bloques básicos. En el presente trabajo se usará para representar el modelo de una motocicleta. Usando este modelo y los resultados de sus simulaciones se recurrirá a una serie de funciones de MATLAB® que serán de gran utilidad a la hora de analizar la controlabilidad del sistema en cuestión. Este es un punto básico de este proyecto, por lo que en esta sección se expondrá el procedimiento que se siguió para la construcción del modelo final incluyendo los ajustes y adaptaciones que se le hicieron al modelo obtenido a partir de las Ecuaciones 15 y 16.

3.1.1. Planteamiento de las ecuaciones para la construcción del modelo

Para empezar la construcción del modelo descrito por la Ecuación 16, lo primero es despejar #� y -�.. Para esto se emplea la Ecuación 20, pues es más sencillo y más organizado dejar estas variables “solas” a un lado de la ecuación bajo la suposición razonable de que la matriz M no es singular.

' #�-�.) = MJKK + MJKB *>?C. ,

Ecuación 20

Para despejar las variables #� y -�. se calcula la inversa de la matriz M, en adelante M-1 (Ecuación 21) por ser la notación estándar.

MJK = LMMMN 1ℎ� + 7ℎ� cos #��ℎOE7��1 − cos� #� + 1 + ℎ� sin #��F 7ℎ� cos #ℎ�E7��1 − cos� #� + 1 + ℎ� sin #��F7ℎ� cos #ℎ�E7��1 − cos� #� + 1 + ℎ� sin #��F 17��1 − cos� #� + 1 + ℎ� sin #�� PQQ

QR

Ecuación 21

Buscando presentar las operaciones y los resultados de la forma más clara y organizada posible se definen cada uno de los elementos de la matriz M-1 en la Ecuaciones 22, 23 y 24. Igualmente se definen los elementos de la matriz K (Ecuación 19) en las Ecuaciones 25 y 26. Esto hará más agradable seguir las operaciones entre matrices.

SJKKK = 1ℎ� + 7ℎ� cos #��ℎOE7��1 − cos� #� + 1 + ℎ� sin #��F

Ecuación 22

16

SJKK� = SJK�K = 7ℎ� cos #ℎ�E7��1 − cos� #� + 1 + ℎ� sin #��F

Ecuación 23

SJK�� = 17��1 − cos� #� + 1 + ℎ� sin #��

Ecuación 24

TKK = �ℎ sin # + 7∆� sin & cos #� + 1 + ℎ� sin #�ℎ�-.� cos #

Ecuación 25

T�K = −2ℎ�-.#�1 + ℎ� sin #� cos # − 7ℎ�#� � sin #

Ecuación 26

Multiplicando la matriz M-1 (Ecuación 21) y la matriz K (Ecuación 19) se obtiene el siguiente resultado:

MJKK= 'SJKKK × TKK + SJKK� × T�KSJK�K × TKK + SJK�� × T�K)

Ecuación 27

Por otro lado el resultado de la multiplicación de la matriz M-1 (Ecuación 21), la matriz B (Ecuación 18) y el vector de entradas del sistema, en ese orden, se presenta en la Ecuación 28.

MJKB *>?C. , = VSJKKK × 7ℎ-. cos # + SJKK� × {−E7�� + ℎ sin #1 + ℎ� sin #�F-.} SJKK� × 1�SJK�K × 7ℎ-. cos # + SJK�� × {−E7�� + ℎ sin #1 + ℎ� sin #�F-.} SJK�� × 1�Y *>?C. ,

Ecuación 28

Reemplazando las Ecuaciones 27 y 28 en la Ecuación 20, se obtienen las igualdades descritas en la Ecuación 29 y en la Ecuación 30, que representan la dinámica del sistema. En estas dos ecuaciones se puede ver la influencia de las entradas >? y C. en el comportamiento dinámico del sistema. A partir de estas dos ecuaciones se empieza a construir el modelo de Simulink™.

#� = Z 1ℎ� + 7ℎ� cos #��ℎOE7��1 − cos� #� + 1 + ℎ� sin #��F[ ��ℎ sin # + 7∆� sin & cos #�+1 + ℎ� sin #�ℎ�-.� cos # �

+ 7ℎ� cos #ℎ�E7��1 − cos� #� + 1 + ℎ� sin #��F <−2ℎ�-.#�1 + ℎ� sin #� cos # − 7ℎ�#� � sin #=

+EZ 1ℎ� + 7ℎ� cos #��ℎOE7��1 − cos� #� + 1 + ℎ� sin #��F[ 7ℎ-. cos #

+ 7ℎ� cos #ℎ�E7��1 − cos� #� + 1 + ℎ� sin #��F −E7�� + ℎ sin #1 + ℎ� sin #�F-.�F>?

+ � 7ℎ� cos #ℎ�E7��1 − cos� #� + 1 + ℎ� sin #��F� 1� C.

Ecuación 29

17

-�. = 7ℎ� cos #ℎ�E7��1 − cos� #� + 1 + ℎ� sin #��F � �ℎ sin # + 7∆� sin & cos #�+1 + ℎ� sin #�ℎ�-.� cos # �

+ 17��1 − cos� #� + 1 + ℎ� sin #�� <−2ℎ�-.#�1 + ℎ� sin #� cos # − 7ℎ�#� � sin #=

+E 7ℎ� cos #ℎ�E7��1 − cos� #� + 1 + ℎ� sin #��F 7ℎ-. cos #

+ 17��1 − cos� #� + 1 + ℎ� sin #�� −E7�� + ℎ sin #1 + ℎ� sin #�F-.�F>?

+ � 17��1 − cos� #� + 1 + ℎ� sin #�� 1� C.

Ecuación 30

Las Ecuaciones 29 y 30 fueron escritas de tal forma que se puedan ver como las Ecuaciones 31 y 32. Al hacer esto se quiso proponer una organización para el modelo de Simulink™ como la que se muestra en la Figura 6. De nuevo, el objetivo es obtener un modelo organizado y fácil de conceptualizar para hacer análisis rápidos y acertados cuando que se requiera.

#� = \K × >? + \� × C. + ]�K

Ecuación 31

-�. = \^ × >? + \O × C. + ]��

Ecuación 32

Figura 6

El modelo de la figura 6 ilustra un sistema en el que H1 representa la función de transferencia entre #� y >? y H3 la función de transferencia entre -�. y >?. Siguiendo ese orden de ideas H2 y H4 son el equivalente entre las salidas y C.. Por su parte NL1 y NL2 representan términos no lineales que se le suman a las salidas y que se eliminarían después de linealizar el sistema (en este caso particular esto no será así, pues todos los términos de la ecuaciones 29 y 30 contienen variables de estado). Esto simplificaría enormemente el problema de control, pues se podrían construir cuatro controladores independientes para manipular la planta.

18

Es importante tener en cuenta que el modelo propuesto en la Figura 6 es un modelo generalizado que en este caso específico no traería las ventajas que se mencionaron, pues los términos que conformarían los bloques NL1 y NL2 no son independientes de las entradas y además contienen variables de estado. Esto hace que el modelo no tenga términos no lineales que desaparezcan después de la linealización. Por otro lado NL1 y NL2 dependen de la entrada >? indirectamente por contener la variable �, que es la integral de >?. Esto implica que la relación entre las entradas y las salidas no se pueden describir únicamente a través de los bloques H1, H2, H3 yH4, ni siquiera después de linealizar el modelo. De cualquier manera se decidió construir el modelo bajo ese formato por considerarlo una forma organizada de distribuir los bloques en Simulink™.

3.1.2. Simplificación del modelo

Antes de construir el modelo se simplificaron las Ecuaciones 29 y 30, que describen la dinámica del sistema, eliminando algunos de sus términos. Para esto se evaluaron expresiones que hacen parte de las ecuaciones mencionadas en sus valores límite (como representación del rango que pueden tomar las variables se incluyó también un valor medio para este análisis) y según el peso en el resultado de cada uno de los términos que las componen, se tomó la decisión de prescindir o no de algunos de ellos. Esto se hizo basándose en la premisa de que entre más simple sea el modelo que se va a trabajar más fácil de entender y de analizar será.

Para darle valor a los parámetros de las ecuaciones se tomó como referencia la moto BMW F 800 GS (ver Anexo 1). Estos mismos parámetros se seguirán usando a lo largo del proyecto. En la Tabla 1 se puede seguir el proceso de simplificación que dio como resultado las Ecuaciones 33 y 34. Como se puede ver, se aprovechó el hecho de que en los valores de límite de la variable #, las funciones sin #� y cos #� son iguales a cero. Esto simplificó bastante el análisis.

� = _ 6a ; ℎ = 0,6; 7 = 0,8; ∆= 0,117; & = 64_ 180a

Ecuación Valores Límite Resultados Numéricos Término que se

conserva

7��1 − cos�#�+ 1 + ℎ� sin #��

# ≈ 0° 0,170� + 1�� = 0 + 1

1 + ℎ� sin #�� # ≈ 45° 0,170,5� + 1,21�� = 0,08 + 1,48

# ≈ 90° 0,171� + 1,31� = 0,17 + 1,73

ℎ sin #+ 7∆� sin & cos #�

# ≈ 0° 0,6 × 0 + 0,04 = 0 + 0,04

ℎ sin # # ≈ 45° 0,6 × 0,7 + 0,03 = 0,42 + 0,03

# ≈ 90° 0,6 × 1 + 0 = 0,6 + 0

−2ℎ�-.#� 1+ ℎ� sin #� cos #− 7ℎ�#� � sin #

# ≈ 0°; #� ≈ i� jk�/m; -. ≈ 66.7�/m −0,62-.#� − 0 × #� � = −64,9 − 0

−2ℎ�-.#� 1+ ℎ� sin #� cos #

# ≈ 45°; #� ≈ i� jk�/m; -. ≈ 66.7�/m −0,54-.#� − 0,17 × #� � = −56,57 − 0,42

# ≈ 90°; #� ≈ i� jk�/m; -. ≈ 66.7�/m 0 × -.#� − 0,25 × #� � = 0 − 0,62

# ≈ 0°; #� ≈ i� jk�/m; -. ≈ 24,17�/m −0,62-.#� − 0 × #� � = −23,54 − 0

# ≈ 45°; #� ≈ i� jk�/m; -. ≈ 24,17�/m −0,54-.#� − 0,17 × #� � = −20,5 − 0,42

# ≈ 90°; #� ≈ i� jk�/m; -. ≈ 24,17�/m 0 × -.#� − 0,25 × #� � = 0 − 0,62

Tabla 1 - Evaluación y eliminación de términos

19

Nótese que para el último término de la Tabla 1 se tuvo en cuenta que -. podría llegar a un máximo de 66,67m/s (240km/h) y que se va a trabajar con un valor promedio de 24,17m/s (87km/h). Adicionalmente

se tiene que, según las simulaciones del modelo que se va a usar (Anexo 2), #� alcanza un valor máximo de i�rad/s cuando # alcanza su valor máximo definido (i�rad). Teniendo en cuenta estos valores en el análisis

se tomó la decisión de conservar únicamente el término que muestra la Tabla 1. A continuación las ecuaciones simplificadas.

#� = Z 1ℎ� + 7ℎ� cos #��ℎO1 + ℎ� sin #��[ �ℎ sin # + 1 + ℎ� sin #�ℎ�-.� cos #�

+ 7ℎ� cos #ℎ�1 + ℎ� sin #�� <−2ℎ�-.#�1 + ℎ� sin #� cos #=

+EZ 1ℎ� + 7ℎ� cos #��ℎO1 + ℎ� sin #��[ 7ℎ-. cos #

+ 7ℎ� cos #ℎ�1 + ℎ� sin #�� −E7�� + ℎ sin #1 + ℎ� sin #�F-.�F>?

+ � 7ℎ� cos #ℎ�1 + ℎ� sin #�� 1� C.

Ecuación 33

-�. = 7ℎ� cos #ℎ�1 + ℎ� sin #�� ℎ sin # + 1 + ℎ� sin #�ℎ�-.� cos #�

+ 11 + ℎ� sin #�� <−2ℎ�-.#�1 + ℎ� sin #� cos #=

+E 7ℎ� cos #ℎ�1 + ℎ� sin #�� 7ℎ-. cos #

+ 11 + ℎ� sin #�� −E7�� + ℎ sin #1 + ℎ� sin #�F-.�F>?

+ � 11 + ℎ� sin #�� 1� C.

Ecuación 34

Complementando lo que se dijo en el primer párrafo de la presente sección, entre más sencillo sea el modelo que refleje el comportamiento del sistema de una forma realista, más fácil será lograr el objetivo de control. Adicionalmente a la hora de construir un controlador real, la simplicidad del modelo que se use puede llegar a afectar positivamente los tiempos de respuesta y por ende el desempeño del controlador.

3.1.3. Adaptación del modelo

Se construyó un modelo en Simulink™ a partir de las ecuaciones originales (Ecuaciones 29 y 30) y otro a partir de las ecuaciones simplificadas (Ecuaciones 33 y 34). Se incluyeron saturadores dentro de cada uno de estos modelos para limitar los valores máximos y mínimos acercando el modelo a la realidad, pues en

20

la práctica las variables tienen limitaciones físicas que se deben tener en cuenta. Se definió un rango de movimiento para � de -30 a 30 grados, un rango de velocidad para -. de 0 a 240km/h y un ángulo de inclinación para # de -90 a 90 grados. Aunque el saturador para # se definió según el rango mencionado anteriormente, para efectos del análisis de controlabilidad se asumirá un rango de -45 a 45 grados, pues se supone que para ángulos mayores la moto caerá al piso. Finalmente se ajustó el modelo de manera que las simulaciones en malla abierta del sistema correspondieran a valores y comportamientos reales de una moto con esas características.

Antes de continuar es importante aclarar que todos los ángulos se trabajaron en radianes por facilidad, pues como se ha mencionado anteriormente, se usó el programa MATLAB® como herramienta a lo largo de este trabajo. Para el resto de medidas se usaron las unidades del SI (Sistema Internacional de Unidades).

Por último se decidió hacer unos cambios en cuanto a las variables de salida. El autor de [5] propone que

se tengan las variables #� y -�. a la salida. Al principio se tomó la decisión de conservar esas dos salidas e

incluir #� , pues se quería tener una mayor información sobre el sistema en este punto. Como variables de

estado se escogieron -., #� y #. De esta forma se empezó a trabajar el sistema, pero llegó un punto (que se discutirá en las conclusiones) en el que se reevaluaron y se cambiaron estas opciones. Finalmente se decidió tener como vector de salida el mismo vector de estado. De esta forma, y como es de esperarse por estar representando un sistema real, se pudo trabajar con una representación en espacio de estados de un sistema estrictamente propio (que su ganancia tiende a cero a medida que se aumenta la frecuencia). Esto se hará evidente más adelante cuando se note la ausencia de la matriz de transmisión directa D en la representación en espacio de estados del sistema.

El modelo correspondiente a las ecuaciones originales se llamó “Modelo_Moto” y el modelo final se llamó “F800GSPS”. Se bautizó con ese nombre por la referencia de la moto (BMW F 800 GS) cuyos parámetros se tomaron para las simulaciones combinado con el hecho de que, por la opción de entradas y salidas que se escogieron, se tuvo que el modelo lineal de la planta G(s) es una matriz rectangular. Esto hace que G(s) no tenga inversa, por lo que se hizo necesario emplear la pseudoinversa para algunos análisis que se llevaron a cabo. Se puede ver el diagrama de bloques con las salidas y las entradas de estos modelos en la Figura 7. Para revisar más a fondo su construcción y sus características se puede ver el Anexo 3.

Figura 7 – Modelos no lineales construidos en Simulink™.

3.1.4. Validación del modelo simplificado

Para validar el modelo se simularon F800GSPS y Modelo_Moto usando las mismas entradas y los mismos parámetros. Para la entrada Fr se usó una señal paso con un valor inicial igual a 0 y un valor final igual a

21

86 que empieza a actuar desde el tiempo t = 0. Para Wsigma se usó una señal seno de amplitud pi/45 (el máximo establecido) y periodo 1/pi que empieza a actuar en el tiempo pi2 segundos. En cuanto a los parámetros, se usaron las especificaciones técnicas de la moto BMW F 800 GS que aparecen en la Tabla 2. En el Anexo 1 se encuentra la hoja de datos de la moto mencionada junto con los cálculos y la explicación de los valores tomados para el modelo. La simulación del modelo correspondiente a las ecuaciones originales (Modelo_Moto) se comparó con la simulación del modelo correspondiente a las ecuaciones simplificadas (F800GSPS). Esto se hizo para validar las simplificaciones que se le hicieron al modelo.

Especificaciones con base en una moto BMW F 800 GS ℎ = 0,6m 7 = 0,8m ∆= 0,117m

& = 64_ 180a rad

� = 300kg

Tabla 2

SIMULACIONES DE Modelo_Moto (izquierda) Y F800GSPS (derecha) PARA LA VALIDACIÓN DE ESTE ÚLTIMO

Fr: Entrada Paso desde t =0s (86N); Wsigma: Entrada seno (f = pi-1rad, Amplitud = pi/45). Empieza a actuar en t = pi2s.

22

Tabla 3

En la Tabla 3 se puede ver que para las mismas entradas (descritas en la parte superior de dicha tabla), la forma de las curvas obtenidas al simular Modelo_Moto son muy similares a las obtenidas al simular F800GSPS. Se puede ver también que las variables de mayor interés, -. y #, alcanzan sus valores máximos en tiempos diferentes para cada uno de los modelos. En las simulaciones de Modelo_Moto se

23

puede ver que # alcanza su valor máximo en el tiempo t = 14.4675s, mientras que en F800GSPS lo alcanza en el tiempo t = 14.4683s. Se presenta una diferencia de 0.8ms. Cabe anotar que el hecho de que # alcance su valor máximo implica que la moto ha caído al piso. Esto hace que la interpretación de los resultados de ahí en adelante sea diferente a la interpretación que se les viene dando hasta ese momento. Por otro lado, en Modelo_Moto, -. alcanza su valor máximo en el tiempo t = 12.7798s, mientras que en F800GSPS lo hace en el tiempo t = 12.7805s. En este caso la diferencia es de 0.7ms. Estas diferencias que se presentan se consideran muy pequeñas si se tiene en cuenta que la moto tardó alrededor de 4.6s en caer al suelo después de haberse presentado una variación en >?.

A pesar de que los resultados obtenidos en las simulaciones de la Tabla 3 no presentan diferencias significativas entre el modelo original y el modelo simplificado, se hicieron otras simulaciones antes de definir qué modelo usar.

En la Tabla 4 se presentan simulaciones de los mismos modelos, pero usando diferentes funciones para las entradas C. y >? (descritas en la parte superior de dicha tabla). De nuevo, se puede ver que las curvas siguen siendo similares. Se presenta una diferencia notoria en -., pero es después de la caída de la moto.

SIMULACIONES DE Modelo_Moto (izquierda) Y F800GSPS (derecha) PARA LA VALIDACIÓN DE ESTE ÚLTIMO Fr: Entrada Paso desde t =0s (86N); Wsigma: Entrada rampa (pendiente = 0,4) hasta t = 10s donde toma un valor constante igual a 4rad/s.

24

Tabla 4

En la Tabla 4 se puede ver que en las simulaciones de Modelo_Moto, # alcanza su valor máximo en el tiempo t = 2.3664s. En las simulaciones de F800GSPS # alcanza su valor máximo en el tiempo t = 2.3843s. Se presenta una diferencia de 17.9ms. Por otro lado, en Modelo_Moto, -. alcanza su valor máximo en el tiempo t = 1.5215s, mientras que en F800GSPS lo hace en el tiempo t = 1.5422s. En este

25

caso la diferencia es de 20.7ms. Aunque las diferencias que se presentan entre una y otra simulación son mayores que en las simulaciones de la Tabla 3, siguen estando en el orden de los milisegundos y se siguen considerando pequeñas al compararlas con los 2.38s que tarda la moto en caerse desde que se presenta un cambio en >?.

En las simulaciones de Modelo_Moto que se muestran en la Tabla 5, # alcanza su valor máximo en el tiempo t = 16.8331s. En las de F800GSPS lo hace en t = 16.8349s. Se presenta una diferencia de 1.8ms. Por otro lado, en Modelo_Moto, -. alcanza su valor máximo en el tiempo t = 14.9035s, mientras que en F800GSPS lo hace en el tiempo t = 14.9047s. Para este caso la diferencia de tiempo que se presenta es de 1.2ms. Al igual que en los casos anteriores, las diferencias que se calculan en este párrafo se pueden considerar muy pequeñas al compararlas con los 6.9s que la moto se demora en caer al piso después de haber modificado el valor de >?.

Hasta ahora se ha podido ver que el sistema simplificado (F800GSPS) tiene tiempos de respuesta más largos que el modelo sin simplificar (Modelo_Moto). Sin embargo al comparar estos retardos con la dinámica de la motocicleta se ha podido establecer que son pequeños. Sólo se presentó una diferencia notoria después de que la moto había caído, por lo que no se tuvo en cuenta para la validación.


Fr: Entrada Paso (86N) desde t = 0s; Wsigma: Entrada paso (0.01rad/s) desde t = pi2s.

26

Tabla 5

Finalmente se presentan las simulaciones de la Tabla 6. En ese caso no se le hacen cambios al valor de la entrada >?. El resultado es igual para los dos modelos. Las curvas tienen la misma forma y las variables toman los mismos valores.

27


Fr: Entrada Paso (86N); Wsigma: Entrada con un valor constante igual a 0rad/s.

Tabla 6

Teniendo en cuenta las pequeñas diferencias que aparecieron entre ambos modelos se toma la decisión de trabajar con el modelo simplificado por considerar que no presenta diferencias significativas con el

28

modelo sin simplificar. Como se dijo antes, esta decisión puede mejorar el tiempo de repuesta de un posible controlador y hacer más sencillo el análisis cuando se están examinando las ecuaciones que describen al sistema. Cada uno de los modelos y su composición detallada se encuentra en el Anexo 3.

En lo que resta de este proyecto se usará el modelo F800GSPS para todos los desarrollos.

3.2. Linealización

En este paso lo primero que se hizo fue escoger el punto de operación alrededor del cual se iba a linealizar el modelo. Teniendo en cuenta que un sistema como el que aquí se viene trabajando puede tener varios puntos de interés alrededor de los cuales linealizar, se evaluaron 21 diferentes. Se debe aclarar que se linealizó alrededor de puntos de operación que se consideraron de interés y que no cumplen con la condición de ser puntos de equilibrio. A estos 21 modelos lineales que se obtuvieron se les hizo el análisis de controlabilidad y posteriormente se seleccionaron algunos para conformar un conjunto de modelos lineales controlables. De lo anterior, se entiende que cada uno de los espacios se genera a partir de diferentes valores de las variables de estado. La Figura 8 pretende ilustrar este punto. Cada una de las imágenes representa un espacio de estados generado a partir de una linealización alrededor de un punto de operación específico. En este ejemplo se quiere mostrar unos espacios de estados controlables (3, 4 y 6) que forman un conjunto de modelos como el que se busca conformar al final del análisis de controlabilidad que se haga en este trabajo.

Figura 8 – Espacios de estado generados a partir de diferentes puntos de operación.

Punto de Operación u0=[Fr,Wsigma] x0=[Vr,theta1,theta]

F800GS_L1 [226,0] [64.2,0,0]

F800GS_L2 [206,0] [58.5,0,0]

F800GS_L3 [186,0] [52.8,0,0]

F800GS_L4 [166,0] [47.15,0,0]

F800GS_L5 [146,0] [41.45,0,0]

F800GS_L6 [126,0] [35.78,0,0]

F800GS_L7 [106,0] [30.1,0,0]

F800GS_L8 [86,0] [24.42,0,0]

F800GS_L9 [66,0] [18.74,0,0]

F800GS_L10 [46,0] [13.06,0,0]

F800GS_L11 [26,0] [7.385,0,0]

F800GS_L12 [6,0] [1.704,0,0]

F800GS_L13 [0,0] [0,0,0]

F800GS_L14 [86,pi/45] [66.67,0.77,0.35]

F800GS_L15 [86,pi/90] [48,0.13,0.065]

29

F800GS_L16 [86,pi/180] [28.1,0.0425,0.028]

F800GS_L17 [86,pi/360] [25.24,0.0192,0.0133]

F800GS_L18 [86,pi/720] [24.62,0.0094,0.0066]

F800GS_L19 [86,pi/1440] [24.48,0.0047,0.0033]

F800GS_L20 [86,pi/2880] [24.44,0.0023,0.0016]

F800GS_L21 [86,pi/5760] [24.42,0.0011,0.0008] Tabla 7 – Puntos de operación alrededor de los cuales se linealizó el modelo.

F800GS_L 17

Se simula la respuesta del sistema F800GSPS a un vector de entrada u0 determinado. Con base en las señales obtenidas en las simulaciones se determina el vector de estados x0 = [Vr,theta1,theta]. Para esto se evalúan las señales que hacen parte de dicho vector en un tiempo t = 12 segundos. Wsigma y Fr son entradas paso. El vector de entrada u0 para esta simulación es :

u0 = [Fr,Wsigma] = [86,pi/360];

RESULTADOS:

La forma de las curvas obtenidas es la esperada. El vector de estado obtenido es el siguiente:

x0 = [25.24,0.0192,0.0133];

Tabla 8 - Simulaciones para la determinación del punto de operación F800GS_L17

Un vector de entradas y un vector de variables de estado constituyen un punto de operación. Una vez más, vale la pena repetir que en este caso se escogieron punto de operación que no cumplen con la condición de ser puntos de equilibrio. Teniendo estos 21 puntos se tomó el modelo F800GSPS y se linealizó alrededor de todos y cada uno de ellos. De aquí en adelante se seguirá y explicará el desarrollo usando como

30

ejemplo los resultados que se obtengan de la linealización del modelo F800GSPS alrededor del punto F800GS_L17. En la Tabla 5 se pueden apreciar las gráficas que se usaron para determinar el punto de operación F800GS_L17. En el Anexo 2 se pueden consultar las otras 20 tablas que contienen la información de los otros 20 puntos.

Para obtener los puntos de operación que se usaron, se simuló la respuesta del modelo a un vector de entrada u0 con valores definidos. Luego se evaluó la respuesta gráficamente para determinar los valores del vector de estados y tener puntos coherentes para la linealización. Para esto se decidió evaluar la respuesta siempre a los 12 segundos de simulación, pues >? empezaba a actuar a los _� segundos y se concluyó que a los 12 segundos había pasado el tiempo suficiente para que el sistema respondiera y a la vez no alcanzara a caer al piso. Siempre se usaron entradas paso para representar a C. y a >?, por lo que su valor a los 12 segundos se conocía de antemano. En la Tabla 8 se pueden ver los resultados de la simulación y el punto de linealización F800GS_L17 resultante.

En el Anexo 4 el lector encontrará el código de MATLAB® empleado para llevar a cabo el procedimiento con cada uno de los 21 puntos mencionados, en todo caso en estas páginas se presentarán todos los resultados relevantes que se obtengan tal y como en la Tabla 7.

D-�.#�#� G = EoF^×^ D-.#�# G + EpF^×� ' C.>?) D�K��^G = EqF^×^ D-.#�# G + ErF^×� ' C.>?)

Ecuación 35

Teniendo los puntos de linealización definidos se usó la función linmod del programa MATLAB® para obtener la linealización del modelo F800GSPS alrededor de los puntos establecidos. Se espera obtener una representación en espacio de estados como la que se muestra en la Ecuación 35. La función linmod arroja como resultado la matriz de estado A, la matriz de entrada B, la matriz de salida C y la matriz de transmisión directa D. Gracias a la escogencia del vector de salidas, que es el mismo vector de estados, se obtuvo una matriz de transmisión directa D igual a cero. Este resultado nos confirma que el sistema que se está representando es estrictamente propio, lo que es de esperarse por tratarse de un sistema real. La representación en espacio de estados del modelo linealizado se llamó F_800GS_L17.

3.3. Controlabilidad

En esta sección se seguirá el proceso para el análisis de controlabilidad de sistemas MIMO propuesto en [1], el libro de control multivariable escrito por Sigurd Skogestad e Ian Postlethwaite, e incluido en el marco teórico del presente trabajo. Se busca determinar si es posible controlar las tres salidas que se tienen. Una vez más, vale la pena aclarar que aunque el procedimiento se hizo para los 21 modelos lineales en las páginas que vienen se explica el procedimiento para uno sólo de ellos. El código de MATLAB® que corresponde a los demás procedimientos se puede consultar en el Anexo 4.

Como anotación final antes de continuar, el análisis en frecuencia para este modelo se hará con un ancho de banda limitado por las características físicas del sistema real. Este ancho de banda está comprendido entre 0rad/s, que corresponde a entradas constantes, y 10rad/s, que es un valor de frecuencia que se ajusta

31

a la realidad; no se espera que bajo condiciones normales un humano acelere y desacelere o gire el timón de un lado al otro más de 1 vez por segundo. Incluso esa frecuencia parece un poco alta.

3.3.1. Escalización

La escalización simplifica el análisis del modelo y obliga al ingeniero a pensar en el desempeño que se espera que tenga el sistema. Para construir las matrices de escalización se debe pensar en las magnitudes aceptables de las variables de entrada, de los disturbios, de los cambios en las referencias y de los errores que se puedan presentar. En resumen, este paso, además de simplificar el diseño de un posible controlador normalizando las variables involucradas en el proceso, ayuda a “aterrizar” el modelo con el que se está trabajando.

rs = 'C.tuv 00 >?tuv) = H226] 00 _45 jk� m⁄ Irx = y-.xtuv 0 00 #�xtuv 00 0 #xtuvz = V66,67 � m⁄ 0 00 _ jk� m⁄ 00 0 _2 jk�Yr. = y-..tuv 0 00 #�.tuv 00 0 #.tuvz = V66,67 � m⁄ 0 00 _ jk� m⁄ 00 0 _4 jk�Y

Ecuación 36

Para escalizar el modelo lineal F_800GS_L17 obtenido en 3.2, lo primero que se hizo fue definir las matrices de escalización, que no son más que matrices diagonales que contienen los valores máximos que se espera tomen las variables de entrada, de salida, los disturbios y las referencias. Como se explica en [1], un importante objetivo del control es minimizar el error {̂ por lo que se escalizó con respecto al error máximo de control. Las matrices de escalización se pueden ver en la Ecuación 40.

De esta forma se obtienen las variables de escalización (Ecuación 41) que se usarán para efectos de control y las funciones de transferencia escalizadas (Ecuación 42) para obtener un modelo en variables escalizadas como el que se presenta en la Ecuación 43. Es importante aclarar que el símbolo ˆ indica que las variables sobre las que aparece no han sido escalizadas.

} = rsJK} ~ , � = rxJK� ~ , { = rxJK{ ~, j = rxJKj ~ Ecuación 37

� = rxJK� � rs

Ecuación 38

� = �}; { = � − j

Ecuación 39

El modelo escalizado se llamó Gs17. Para obtener este resultado se usó el programa MATLAB®.

32

3.3.2. Realización Mínima

Después de escalizar el modelo se debe obtener la realización mínima del mismo. La realización mínima es la representación en espacio de estados que describe el comportamiento de entrada-salida del sistema con el menor número de estados posible. Tener la representación mínima garantiza que se está trabajando con un modelo observable y controlable.

En este caso se pudo comprobar que el modelo Gs17 obtenido en 3.3.1. corresponde a la realización mínima, por lo que se puede seguir trabajando con el resultado obtenido en esa sección. Sin embargo, en las siguientes secciones se decidió trabajar con el resultado (minGs17) obtenido en 3.2.2. para evitar confusiones, pues se pueden presentar casos en los que el resultado obtenido en este punto sea distinto del que se obtuvo en el punto anterior. En la Figura 9 se presenta la representación en espacio de estados de minGs17 dentro de una estructura sys en MATLAB® obtenida usando el código que se puede consultar en el Anexo 4. La columna del extremo derecho y la última fila sirven únicamente para especificar que se trata de una estructura sys de MATLAB®, no hacen parte de la representación en espacio de estados del sistema.

Figura 9

En cuanto a las representaciones en espacio de estado de los modelos linealizados alrededor de los otros puntos se encontró que la realización mínima se había encontrado en el paso anterior al igual que en este caso. Sólo en el modelo que corresponde al punto de linealización F800GS_L13 se obtuvo una respuesta distinta. Esto se explica porque en el punto de linealización F800GS_L13 se tiene el vector de entradas en cero, al igual que el de variables de estado. Por la condición no holonómica del sistema este caso no tiene sentido. Únicamente se usará como ejemplo a lo largo de este trabajo. En la Figura 10 se muestra que para obtener la realización mínima de este sistema fue necesario remover tres estados del mismo. Adicionalmente las únicas matrices cuyos valores son diferentes de cero son B y C, pero sus dimensiones cambian. Esto quiere decir que el modelo linealizado alrededor del punto F800GS_L13 no es controlable, pues no es posible modificar sus estados.

Como se mencionó en un párrafo anterior, de aquí en adelante se trabajará con la representación en espacio de estados minGs17, que es el resultado obtenido en el presente numeral.

33

Figura 10

3.3.3. Controlabilidad funcional

DIRECCIÓN DE SALIDAS EN LA QUE NO HAY GANANCIA - F800GS_L17

DIRECCIÓN DE SALIDA NO CONTROLABLE DIRECCIÓN 1 – -.

DIRECCIÓN 2 – #� DIRECCIÓN 3 – #

Tabla 9

-4

-2

0

2

4

x 10-17

-0.015

-0.01

-0.005

0

0

0.2

0.4

0.6

0.8

1

0 1 2 3 4 5 6 7 8 9 10

-2.5

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

x 10-17 Dirección sin ganancia U17(1,3)

rad/s

U17(1

,3)

0 1 2 3 4 5 6 7 8 9 10-0.012

-0.01

-0.008

-0.006

-0.004

-0.002

0Dirección sin ganancia U17(2,3)

rad/s

U17(2

,3)

0 1 2 3 4 5 6 7 8 9 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


rad/s

U(3

,3)

34

La controlabilidad funcional es una medida de la habilidad de las entradas para afectar las salidas de forma independiente. De acuerdo con lo definido por los autores Sigurd Skogestad e Ian Postlethwaite en su libro “Multivariable Feedback Control: Analysis and Design”, se podría decir que este sistema no es controlable funcionalmente con sólo inspeccionar sus entradas y sus salidas, pues tiene menos entradas que salidas. Esto quiere decir que habrá una dirección de salida en la que las entradas no podrán afectar a las salidas. Esa dirección varía con la frecuencia y se puede ver graficando, en este caso, la tercera columna de la matriz de salidas U que se obtiene al hacer la descomposición en valores singulares del sistema. Esto es evidente si se tiene en cuenta que para un sistema de estas características la tercera fila de la matriz de valores singulares es cero es decir que no hay ganancias en esa fila.

DIRECCIÓN DE SALIDAS EN LA QUE NO HAY GANANCIA - F800GS_L4

DIRECCIÓN DE SALIDA NO CONTROLABLE DIRECCIÓN 1 - -.

DIRECCIÓN 2 - #� DIRECCIÓN 3 - #

Tabla 10

El código de MATLAB® que se implementó para obtener y graficar las direcciones en las que no hay ganancia (Tabla 9), se puede consultar en el Anexo 4.

En la Tabla 9 se pueden ver las direcciones de las salidas del modelo minGs17 (obtenido al linealizar alrededor del punto F800GS_17 y obtener la realización mínima del modelo resultante) en las que no hay ganancia, es decir las direcciones en las que no se puede ejercer control sobre las salidas. Se ve también cada una de sus componentes en frecuencia. Como es de esperarse, sobre la salida -. no se puede ejercer

-1

-0.5

0

0.5

1

-6

-4

-2

0

x 10-3

0

0.2

0.4

0.6

0.8

1

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

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8


rad/s

U4(1

,3)

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

-5

-4

-3

-2

-1

0x 10

-3 Dirección sin ganancia U4(2,3)

rad/s

U4(2

,3)

0 1 2 3 4 5 6 7 8 9 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


rad/s

U4(3

,3)

35

control mientras esté sobre el cero. Esto no es un inconveniente, pues cero no hace parte de los valores aceptados de -.. Lo anterior se concluye desde la sección 2.3 (restricción 1) cuando se especifica que el sistema es no holonómico (si la velocidad de la motocicleta es cero no se puede modificar su orientación y por consiguiente no se puede ejercer control sobre el sistema) y de los resultados que se muestran en la Figura 10 cuando se intenta obtener la realización mínima del modelo linealizado alrededor de el punto F800GS_13 y se llega a que no es posible obtener información sobre las variables de estado y tampoco se pueden modificar (no es observable ni controlable) cuando la velocidad longitudinal es cero.

Para el ángulo de inclinación #, la dirección en la que no hay ganancia tiende a acercarse a 1, que corresponde a un valor de 57.3°, a medida que aumenta la frecuencia. Se espera no llegar a este valor de #, pues se considera que para valores mayores a 45° la moto caerá al suelo irremediablemente. Los demás valores de # en los que no hay ganancia se deben evitar. Esto no representa un problema porque no se espera que # mantenga una oscilación a una frecuencia constante con una amplitud determinada durante mucho tiempo y esto implica que los valores en los que no se puede controlar sean pasajeros. Adicionalmente se ve que en DC no se puede controlar # cuando es cero. Esto último tampoco representa un problema, pues mientras # sea cero, el objetivo de control, que es evitar que la motocicleta caiga al piso, se está cumpliendo.

En cuanto a la velocidad de inclinación #� se puede decir, con base en su gráfica, que los valores en los que no se puede manipular son cercanos a cero sin llegar a este valor, al menos en las frecuencias de interés de

este sistema. Por ser un sistema tan inestable no es factible que #� tome alguno de esos valores durante un

tiempo suficiente para afectar la controlabilidad de esta salida. Los valores de #� que se muestran en la gráfica se deben evitar.

La gráfica tridimensional muestra las direcciones de salida que se deben evitar para que las variables de salida siempre sean controlables.

La Tabla 10 contiene la información de las direcciones de las salidas del modelo minGs4 (obtenido al linealizar alrededor del punto F800GS_4) en las que no hay ganancia. La forma de las curvas es la misma para los dos modelos, e incluso se llegan a las mismas conclusiones al observar las gráficas.

3.3.4. Obtención de los polos del sistema

Para encontrar los polos del sistema se usó la función opde creada por Kjetil Havre para [1], el libro de control multivariable escrito por Sigurd Skogestad e Ian Postlethwaite. Esta función encuentra los polos de un sistema y sus direcciones. El código de MATLAB® que se usó para esto se puede ver en el Anexo 4.

En la Tabla 11 se presentan los polos de cada uno de los modelos y sus direcciones. En la primera columna de esta tabla aparece el nombre del modelo al que corresponden los polos obtenidos, en la segunda aparecen los polos del sistema organizados en un vector de dimensiones 4x1 y en la última se muestran las direcciones de cada polo organizadas en una matriz de dimensiones 3x4. La primera columna de la matriz de direcciones es un vector que corresponde a la dirección del polo de la primera fila del vector de polos que se muestra en la segunda columna de la Tabla 11. La segunda columna de la matriz de direcciones corresponde a la dirección del polo de la segunda fila del vector de polos y así sucesivamente. Como se puede ver la matriz de direcciones tiene cuatro columnas en las que se encuentran las direcciones de los cuatro polos que se obtuvieron.

36

Se puede ver que los polos de los primeros 13 modelos son iguales tanto en magnitud como en dirección. Del modelo número 14 en adelante se empiezan a presentar diferencias, sin embargo del modelo 16 al 21 los polos -0.13 y 0.13 se mantienen, pero con cambios en su dirección con respecto a los anteriores. Estos resultados indican que si >? se mantiene en cero, C. puede tomar diferentes valores dentro del punto de operación alrededor del cual se linealiza sin que ello genere un vresultado diferente en cuanto a los polos se refiere.

MODELO POLOS DIRECCIONES DE LOS POLOS

minGs1 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

minGs2 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

minGs3 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

minGs4 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

minGs5 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

minGs6 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

minGs7 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

minGs8 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

minGs9 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

minGs10 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

minGs11 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

minGs12 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

minGs13 y−0.13000.13 z D 0 1−0.06 00.99 0 0 00 0.06−1 0.99G

37

minGs14 y−0.1800.030.05 z D−0.98 −0.970.02 0−0.18 0.22 −0.96 −0.944 × 10J^ 0.010.29 0.33 G

minGs15 y−0.13−0.0600.1 z D 0.95 0.97−0.02 6.3 × 10J^0.28 −0.19 0.82 −0.560 0.04−0.56 0.82 G

minGs16 y −0.1308.9 × 10JO0.13 z D−0.36 −0.990.06 0−0.93 0.04 0.99 −0.357−2.6 × 10J� 0.06−0.06 0.93 G

minGs17 y −0.1301.7 × 10JO0.13 z D−0.17 −0.990.06 0−0.98 0.02 0.99 −0.17−2.6 × 10J� 0.06−0.03 0.98 G

minGs18 y −0.1304 × 10J�0.13 z D−0.08 −0.990.06 0−0.99 0.01 0.99 −0.08−3 × 10J� 0.06−0.02 0.99 G

minGs19 y −0.13−3 × 10J^00.12 z D−0.11 0.990.06 1.1 × 10J�−0.99 −7.5 × 10J^ −0.91 −0.10 0.06−0.41 0.99G

minGs20 y −0.1302.7 × 10J�0.13 z D−0.02 −10.06 0−0.99 0.02 1 −0.02−5 × 10J� 0.06−3.7 × 10J^ 0.99 G

minGs21 y −0.1306.6 × 10J�0.13 z D−0.01 −10.06 0−0.99 1.3 × 10J^ 1 −0.01−6.2 × 10JK� 0.06−1.8 × 10J^ 0.99 G

Tabla 11 Lista de los polos de todos los modelos y sus direcciones

También se puede ver en la Tabla 11 que no se obtuvieron polos complejos y que los polos que aparecen en el semiplano derecho son cercanos al origen, por lo que a pesar de ser inestables no representan un problema serio a la hora de implementar una acción de control. Adicionalmente en los primeros 13 modelos hay polos negativos y positivos alrededor del eje imaginario y al observar sus direcciones se puede establecer que se presenta una simetría con respecto a un plano.

3.3.5. Obtención de los ceros del sistema

Para encontrar los ceros del sistema se usó la función ozde creada por Kjetil Havre para [1], el libro de control multivariable escrito por Sigurd Skogestad e Ian Postlethwaite. Esta función encuentra los ceros de un sistema y sus direcciones. El código de MATLAB® que se usó para esto se puede ver en el Anexo 4.

MODELO CEROS

minGs1 -80.26

minGs2 -73.13

minGs3 -66.01

minGs4 -58.94

minGs5 -51.82

minGs6 -44.73

minGs7 -37.63

minGs8 -30.53

38

minGs9 -23.43

minGs10 -16.33

minGs11 -9.24

minGs12 -2.14

minGs13 -

minGs14 -83.32

minGs15 -59.99

minGs16 -35.13

minGs17 -31.56

minGs18 -30.78

minGs19 -83.33

minGs20 -30.56

minGs21 -30.53

Tabla 12 Lista de los ceros de todos los modelos

Como se puede ver en la Tabla 10 todos los ceros de todos los modelos (esto no se cumple para el F800GS_L13 que es un caso especial que no es coherente con las restricciones del sistema) están en el semiplano izquierdo, por lo que en este sentido no se presenta dificultad alguna y no es necesario conocer sus direcciones.

3.3.6. Matriz RGA y su respuesta en frecuencia

En la Figura 12 se presentan los valores en frecuencia de la RGA correspondientes al modelo minGs17. Como se ve en estas gráficas, los valores de la RGA nunca son mayores a 1. Esto es una buena noticia si se quiere implementar un controlador para este sistema en el futuro. Las gráficas que se encuentran en la tabla de la Figura 12 contienen información sobre la relación que existe entre las entradas (representadas por las columnas de la tabla) y las salidas (representadas por las filas de la tabla). En la gráfica “RGA (1,1) vs w” se puede ver que para todos los valores de frecuencia, la salida -. depende exclusivamente de la entrada C.. Por otro lado en la gráfica “RGA (3,2) vs w” se ve que a bajas frecuencias, que son las de interés para este trabajo, la salida # depende en gran medida de la entrada >?, mientras que en la gráfica

“RGA (2,2) vs w” se ve que a estas frecuencias la influencia de le entrada >? sobre la salida #� es muy débil. Al aumentar la frecuencia esta situación se revierte al punto que >? deja de tener influencia sobre #

y presenta una influencia muy fuerte sobre #� . Es importante señalar que a frecuencias altas ninguna de las entradas tiene influencia sobre # y por consiguiente no es posible controlar esa salida.

Figura 11

39

Se quiso evaluar la RGA en DC por lo que se tomó la decisión de evaluarla en una frecuencia igual a1e-6 rad/s, que es muy cercano a cero y se encuentra dentro del ancho de banda. El resultado es la matriz RGA que se presenta en la Figura 11. De acuerdo con esta matriz, a muy bajas frecuencias, la velocidad lineal

de la moto se controlará exclusivamente con la entrada C., la velocidad angular #� será la más difícil de controlar y se deberá hacer mediante la entrada >?. Esta última entrada también será la que controle a #, que será relativamente fácil de controlar.

C. >?

-.

#�

#

Figura 12 - Valores de la RGA en función de la frecuencia w (rad/s)

Es importante resaltar que la respuesta en frecuencia de los valores de la RGA es coherente con los resultados obtenidos en la sección 3.3.5., donde se determinó que todos los ceros del sistema se encontraban en el semiplano izquierdo. Lo anterior, porque cuando algún elemento de la RGA cambia de signo entre 0 e infinito, quiere decir que hay al menos un cero en el semiplano derecho.

0 1 2 3 4 5 6 7 8 9 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1RGA(1,1) vs w

0 1 2 3 4 5 6 7 8 9 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1RGA(1,2) vs w

0 1 2 3 4 5 6 7 8 9 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1RGA(2,1) vs w

0 1 2 3 4 5 6 7 8 9 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1RGA(2,2) vs w

0 1 2 3 4 5 6 7 8 9 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1RGA(3,1) vs w

0 1 2 3 4 5 6 7 8 9 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1RGA(3,2) vs w

40

Viendo la matriz de la RGA se puede determinar si se está trabajando con entradas de más. En este caso los elementos de las columnas suman 1. Esto quiere decir que el número de entradas es suficiente y que no se tienen entradas de más. Si los elementos de las columnas suman mucho menos que 1, quiere decir que se puede prescindir de una o más entradas que no tienen mucha influencia sobre las salidas. Por otro lado cuando los elementos de las filas no suman 1 quiere decir que hay salidas que no se pueden controlar. En este caso, a bajas frecuencias dos de las salidas (-. y #) se pueden controlar, una mejor que la otra, y la

salida restante (#�) presenta problemas. A altas frecuencias las dos primeras (-. y #�) se pueden controlar y la restante (#) no. Eso viendo únicamente la RGA, pero hay que aclarar que al cruzar esta información con los valores singulares se ve que hay frecuencias altas a las que es muy difícil ejercer control.

3.3.7. SVD y respuesta en frecuencia de los valores singulares

En la Figura 13 se puede ver la respuesta en frecuencia de los valores singulares del modelo minGs17. Las frecuencias de interés de este sistema son bajas, por esta razón las simulaciones que se han hecho van hasta una frecuencia máxima de 10rad/s. Observando la simulación de los valores singulares se puede ver que a partir de los 0.08rad/s, el valor singular mínimo �� correspondiente a >? (en rojo) se vuelve

menor a 1 y a partir de los 0,34rad/s pasa lo mismo con el valor singular máximo �� que le corresponde a C. (en amarillo). A frecuencias en las que los valores singulares son menores a 1 se puede presentar saturación, pues será necesario usar entradas muy grandes para obtener los cambios deseados a la salida. Como minGs17 es un modelo escalizado, se podrá lograr una variación mínima de �� a la salida en

cualquier dirección con una entrada de valor unitario. Por esta razón a frecuencias a las que se quiere ejercer control se busca que como mínimo �� sea mayor a la unidad para evitar que el sistema se sature.

Figura 13 – Respuesta en frecuencia de los valores singulares de minGs17

Por otro lado el número de condición ��, cuya definición se muestra en la Ecuación 44, toma valores grandes (mayores a 10) en este rango de frecuencias. Por lo general esto se asocia con una alta sensibilidad a las incertidumbres y se dice que la planta está mal condicionada, sin embargo esto no siempre es cierto ya que se pueden presentar casos, como el que se está analizando, en los que el número

41

de condición sea grande porque el valor singular máximo es muy grande. Mientras el valor singular mínimo sea mayor a 1, tener un número de condición grande no representará un problema. En este caso el valor singular mínimo es mayor a 1 para frecuencias menores o iguales a 0,076rad/s, para frecuencias mayores se pueden presentar problemas de controlabilidad.

�� = �� ⁄

Ecuación 40

42

4. CONTROLADOR PROPUESTO

Para el sistema minGs17 se necesita un controlador que mantenga el ángulo # en cero (o un valor deseado cercano a cero) y la velocidad lineal en el valor que se desee (dentro del rango establecido para esta variable). En esta sección se tendrán en cuenta los resultados del análisis de controlabilidad y se propondrá, si de acuerdo con la información que se obtuvo del análisis es posible, un controlador sencillo.

Tanto la Figura 11 como la Figura 12 muestran claramente que la salida -. del sistema minGs17 se puede controlar usando únicamente la entrada C.. Por esta razón se ha decidido diseñar un controlador SISO para controlar la salida mencionada.

Usando MATLAB® se obtuvo la función de transferencia Gfv que describe la relación entre las variables C. y -.. Esta ecuación se factorizó como se ve en la Ecuación 41 para facilitar su análisis.

��- = 80.23 × 10J^ m127.84 × 10J^ − m�127.84 × 10J^ + m�m128.18 × 10J^ − m�127.51 × 10J^ + m�173.43 × 10J� + m�

Ecuación 41

Al examinar la Ecuación 41, se encontraron términos en el numerador muy similares a algunos del denominador. Se decidió cancelar entre sí estos términos y se obtuvo lo siguiente:

��- = 80.23 × 10J^ 1173.43 × 10J� + m

Ecuación 42

Al ver la Ecuación 42 se decidió diseñar el controlador SISO usando el método analítico de síntesis directa para sistemas de primer orden. Para esto se reescribe la Ecuación 42 para organizar los términos como se muestran a continuación:

�� = �/�m + 1 .x�.�u��u�� 6�� - = 462.65766m + 1

Ecuación 43

Teniendo la función de transferencia escrita como en la Ecuación 43 se diseñó un controlador PI por síntesis directa usando las definiciones que se conocen y los valores específicos del sistema que se viene trabajando (ver Ecuación 44).

�9m� = �9 �1 + 1��m� = ��/� ��m + 1�m � = 5766�462.6 �5766m + 15766m � = 100 �5766m + 15766m �

Ecuación 44

En la Ecuación 44 se puede cambiar la velocidad de respuesta del sistema variando �, pues esta variable corresponde a la velocidad de respuesta del mismo. En este caso se escogió � = 0.1246, aproximadamente. En la Figura 14 se puede ver la respuesta usando el controlador diseñado con el sistema descrito por la función de transferencia que aparece en la Ecuación 41. En esa figura se muestra cómo la salida -. del modelo minGs17 es controlada utilizando únicamente la entrada C..

43

Figura 14

La respuesta de este sistema presenta un sobrepico del 0.001%, llega al valor de la referencia (40) en 0.865s y tiene un error en estado estacionario del 0.00016% (ver Figura 15). Este es un desempeño aceptable para la velocidad de una motocicleta.

44

Figura 15 – Respuesta de vr cuando se implement el controlador propuesto.

Para las otras dos salidas no se obtuvo un resultado como el que se obtuvo para -.. Los resultados de la RGA no indican una relación entre estas salidas y alguna de las entradas que permita concluir que un controlador SISO sea la solución adecuada para su control.

45

5. ANÁLISIS DE RESULTADOS Y CONCLUSIONES

Para el modelo minGs17 con el que se viene trabajando se puede ver que en las frecuencias más bajas los valores singulares (Figura 16) son mayores a 1 e incluso el valor singular mayor alcanza una ganancia de 664 en 0.01rad/s, como se aprecia en la figura anterior. En ese mismo punto el valor singular menor tiene una ganancia de 3.467. Esto hace que en este punto el número de condición sea grande (191,52), sin embargo, al encontrarse los dos valores singulares por encima de 1 no se presentarán problemas de saturación por entradas muy grandes.

Figura 16

El valor singular menor alcanza su frecuencia de corte con 1 alrededor de los 0,076rad/s. En este punto el valor del número de condición es 28,44. En el valor de corte con 1 del primer valor singular el número de condición es 4,47. Finalmente en 10rad/s, el número de condición tiene un valor de 31,79. En este rango de frecuencias el número de condición es significativamente mayor en las frecuencias más bajas, sin embargo mientras los valores singulares sean mayores a 1 no se presentarán problemas relacionados con este resultado. Adicionalmente en las frecuencias más altas, donde el número de condición no es tan

46

grande, sí se presenta un problema de controlabilidad, pues la ganancia de los valores singulares es demasiado pequeña y esto hace que la influencia de las entradas sobre las salidas sea muy débil.

Como se ve en la Figura 12, los valores de la RGA nunca son mayores a 1 en ningún valor de frecuencia. Las gráficas que se encuentran en la tabla de la Figura 12 contienen información sobre la relación que existe entre las entradas (representadas por las columnas de la tabla) y las salidas (representadas por las filas de la tabla). En la gráfica “RGA (1,1) vs w” se puede ver que para todos los valores de frecuencia, la salida -. depende exclusivamente de la entrada C.. Por otro lado en la gráfica “RGA (3,2) vs w” se ve que a bajas frecuencias, que son las de interés para este trabajo, la salida # depende en gran medida de la entrada >?, mientras que en la gráfica “RGA (2,2) vs w” se ve que a estas frecuencias la influencia de le

entrada >? sobre la salida #� es muy débil. Al aumentar la frecuencia esta situación se revierte al punto

que >? deja de tener influencia sobre # y presenta una influencia muy fuerte sobre #� . Es importante señalar que a frecuencias altas ninguna de las entradas tiene influencia sobre # y por consiguiente no es posible controlar esa salida a esas frecuencias.

Viendo la matriz de la RGA a bajas frecuencias (Figura 11) se puede determinar si se está trabajando con entradas de más. En este caso los elementos de las columnas suman 1. Esto quiere decir que el número de entradas es suficiente y que no se tienen entradas de más. Por otro lado cuando los elementos de las filas no suman 1 quiere decir que hay salidas que no se pueden controlar. En este caso, a bajas frecuencias dos

de las salidas (-. y #) se pueden controlar, una mejor que la otra, y la salida restante (#�) presenta

problemas. A altas frecuencias las dos primeras (-. y #�) se pueden controlar y la restante (#) no. Lo anterior se concluye de observar únicamente la RGA, pero también se deben tener en cuenta los demás resultados que se obtuvieron a lo largo del análisis. En este caso particular, la gráfica de los valores singulares muestra que a frecuencias altas es muy difícil ejercer control debido a la alta sensibilidad a las incertidumbres que presenta el sistema y a los posibles problemas de saturación que se tendrán por la magnitud de las entradas necesarias para ejercer control sobre las salidas.

De acuerdo con la matriz de la Figura 11, a muy bajas frecuencias, la velocidad lineal -. de la moto se

controlará exclusivamente con la entrada C. y la velocidad angular #� será la más difícil de controlar y se deberá hacer mediante la entrada >?. Esta última entrada también servirá para controlar a #, que será relativamente fácil de controlar. Sin embargo, como se advirtió en la sección 4, esta última relación no es

clara, por lo que no se puede afirmar que las salidas # y #� dependan exclusivamente de la entrada >? en ningún caso.

Algunos sistemas MIMO se pueden controlar con controladores SISO. En este punto el análisis que se haga soportado por la matriz RGA será de gran ayuda. El análisis en frecuencia de sus elementos entrega información sobre la relación que existe entre las entradas y las salidas. En este caso específico, se pudo concluir que a bajas frecuencias, que en realidad son las de interés en este trabajo, se puede controlar una de las salidas del sistema (-.) con un controlador SISO. A altas frecuencias el sistema es muy difícil de controlar, sin embargo esto no representa un problema porque no se espera tener entradas a estas frecuencias, sólo en eventualidades como cuando se golpea una piedra. De cualquier forma existen sistemas de amortiguación para el timón diseñados para que esta situación no se presente.

El hecho de que el sistema no sea controlable funcionalmente no quiere decir que no se pueda controlar, pero sí que no se pueden controlar las salidas independientemente. Esto hará que, por ejemplo, al cambiar

el valor de # se presenten variaciones en el valor de #� . Lo anterior es lógico por la relación diferencial que existe entre las dos variables.

47

El conjunto de modelos controlables del que se habló está conformado por todos los modelos exceptuando el que se obtuvo a partir del punto de linealización F800GS_L13. Esto es de esperarse, pues se sabe que una moto detenida es muy difícil de controlar.

Se sintetizó un controlador PI de una entrada y una salida por síntesis directa para controlar -.. EL desempeño del controlador propuesto es aceptable, de acuerdo con lo manifestado en la sección 4. Para las otras dos salidas no se obtuvo un resultado como el que se obtuvo para -.. Los resultados de la RGA no indican una relación entre estas salidas y alguna de las entradas que permita concluir que un controlador SISO sea la solución adecuada para su control, por lo que probablemente sea necesario implementar una estrategia de control más complicada que no hace parte del alcance del presente trabajo.

48

6. BIBLIOGRAFÍA

[1] SKOGESTAD, Sigurd y POSTLETHWAITE, Ian. Multivariable Feedback Control: Analysis and Design. 2 ed. Chichester: John Wiley & Sons, 2005. 608p.

[2] LIMEBEER David J.N. y SHARP, Robin S. Single-Track Vehicle Modeling and Control: Bicycles, Motorcycles, and Models. En: IEEE Control Systems Magazine. oct. 2006; p. 34 – 61.

[3] YI, Jingang et al. Trajectory Tracking and Balance Stabilization Control of Autonomous Motorcycles. En: INTERNATIONAL CONFERENCE ON ROBOTICS AND AUTOMATION (5 : 2006 : Orlando, Florida). Proceedings of the 2006 IEEE. Orlando: IEEE, 2006.

[4] SPONG, Mark W. et al. Robot Modeling and Control. 1ed. Chichester: John Wiley & Sons, Inc., 2005. 407p.

[5] GETZ, Neil H. Dynamic Inversion of Nonlinear Maps with Applications to Nonlinear Control and Robotics. A dissertation submitted in partial satisfaction of the requirements for the degree of Doctor of Philosophy in Electrical Engineering and Computer Sciences in the GRADUATE DIVISION of the UNIVERSITY of CALIFORNIA at BERKELEY. 1995.

[6] http://www.mathworks.com/access/helpdesk/help/toolbox/simulink/slref/linmod2.html

[7] http://www.mathworks.de/access/helpdesk_r13/help/toolbox/mutools/pck.html

[8] http://www.mathworks.com/access/helpdesk_r13/help/toolbox/mutools/veval.html

[9] http://www.in.gov/history/2918.htm

[10] http://www.ozebook.com/roper.htm

[11] http://www.vf750fd.com/blurbs/first.html

[12] COSSALTER, Vittore. Motorcycle Dynamics. Greendale, WI: Race Dynamics, 2002.

49

7. ANEXOS

50

ANEXO 1

51

Para el desarrollo del presente trabajo se tomó como referencia la moto trail BMW F 800 GS. Se puede decir que las motos trail son una versión más versátil de las motos de enduro o de cross, pues además de tener un buen desempeño en terrenos difíciles cuentan con el equipamiento adecuado para viajar con seguridad por una carretera. A continuación se muestran unos recortes tomados de la hoja de datos publicada por BMW en su página web http://www.bmwmotorcycles.com/us/en/index.html .

De esta información se tomó la distancia ∆ entre el punto de contacto de la rueda delantera y la proyección del tenedor delantero sobre el suelo, que aparece en esta hoja de datos como “Castor” con un valor de 117mm. También se tomó el ángulo η del tenedor delantero con respecto a la horizontal, cuyo valor es de 64° y aparece bajo el nombre de “steering head angle”. En cuanto a la masa, se tomó la decisión de usar 300kg, ya que es una opción razonable dentro de un rango que va desde 207kg (con el tanque lleno, pero sin ningún otro peso adicional) y 443kg (que es el peso máximo permitido).

Los parámetros restantes no aparecen explícitamente dentro de estos datos, razón por la cual se calcularon o se supusieron con base en datos que sí aparecen en este cuadro. El primero es el caso de la fuerza de tracción longitudinal de la llanta trasera Fr y el segundo es el de las coordenadas del centro de masa b

52

(Distancia horizontal entre el punto de contacto de la rueda trasera y la proyección del centro de masa sobre el suelo) y h (la altura del centro de masa de la motocicleta).

Para calcular Fr se supuso que el torque máximo (83N) que aparece en esta hoja de datos se ejerce sobre la llanta trasera. Se tiene además que el rin tiene un diámetro de 17” que equivale a 0,4318m. Se tiene también que, de acuerdo sus especificaciones, el neumático trasero tiene una sección transversal de 150mm. Si se supone que es cilíndrico y que la parte que queda cubierta por el rin es despreciable se puede decir que el radio de la llanta, incluyendo el neumático y el rin, es el resultado de la Ecuación A1 – 1.

j = 0,4318� + 0,15 × 2��2 = 0,3659�

Ecuación A1 - 1

Por otro lado se sabe que el torque se define como

� = j × C

Ecuación A1 - 2

Reemplazando por los valores que se tienen y despejando se llega a que C. = 226,83]. En el proyecto se tomará como 226N, sin decimales.

Finalmente para el centro de masa se tuvo en cuenta que la distancia entre las dos llantas o wheelbase, como aparece en la tabla, es de 1,578m. Se decidió escoger b como la mitad (aproximadamente) de esta longitud y se determinó que el valor de b sería 0,8. Para h se decidió que sería igual al 70% de la altura de la silla más baja (850mm) por lo que se determinó que el valor de h sería 0.6

A continuación se presenta el resumen de los valores obtenidos a partir de la hoja de datos de la moto BMW F 800 GS en la Tabla A1 - 1. Al igual que en el resto del presente trabajo, los valores de los parámetros que aparecen en la tabla están expresados en unidades del SI (Le Système International d'Unités), excepto por los grados que se expresarán en radianes.

Especificaciones con base en una moto BMW F 800 GS ℎ = 0,6m

7 = 0,8m

∆= 0,117m

& = 64_ 180a rad

� = 300kg

Tabla A1 - 1

53

ANEXO 2

F800GS_L1

Se simula la respuesta del sistema “F800GSPS” a un vector de entrada u0 determinado. Con base en las señales obtenidas en las simulaciones se determina el vector de estados x0. Para esto se evalúan las señales que hacen parte de dicho vector en un tiempo t = 12 segundos. El vector de entrada u0 para esta simulación, en t = 12 segundos, es :

u0=[Fr,Wsigma]=[226,0];

RESULTADOS:

Observando la simulación se puede decir que con las entradas que se usaron la moto anda en línea recta sin caer al piso durante los 30 segundos que dura la misma y que su velocidad aumenta de forma lineal hasta llegar a la velocidad máxima que se determinó cuando se adaptó el modelo a los requerimientos del proyecto.

El vector de estados obtenido a los 12 segundos de simulación es el siguiente:

x0=[Vr,theta1,theta]=[64.2,0,0];

54

F800GS_L2


u0=[Fr,Wsigma]=[206,0];

RESULTADOS:




55

F800GS_L 3


u0=[Fr,Wsigma]=[186,0];

RESULTADOS:




56

F800GS_L 4


u0=[Fr,Wsigma]=[166,0];

RESULTADOS:



x0=[Vr,theta1,theta]=[47.15,0,0];

57

F800GS_L 5


u0=[Fr,Wsigma]=[146,0];

RESULTADOS:



x0=[Vr,theta1,theta]=[41.45,0,0];

58

F800GS_L 6


u0=[Fr,Wsigma]=[126,0];

RESULTADOS:



x0=[Vr,theta1,theta]=[35.78,0,0];

59

F800GS_L 7


u0=[Fr,Wsigma]=[106,0];

RESULTADOS:




60

F800GS_L 8


u0=[Fr,Wsigma]=[86,0];

RESULTADOS:

Observando la simulación se puede decir que con las entradas que se usaron la moto anda en línea recta sin caer al piso durante los 30 segundos que dura la misma y que su velocidad aumenta de forma lineal con el tiempo. Con el torque aplicado no alcanza la velocidad máxima en el tiempo que dura la simulación.


x0=[Vr,theta1,theta]=[24.42,0,0];

61

F800GS_L 9



RESULTADOS:



x0=[Vr,theta1,theta]=[18.74,0,0];

62

F800GS_L 10



RESULTADOS:



x0=[Vr,theta1,theta]=[13.06,0,0];

63

F800GS_L 11



RESULTADOS:



x0=[Vr,theta1,theta]=[7.385,0,0];

64

F800GS_L 12



RESULTADOS:



x0=[Vr,theta1,theta]=[1.704,0,0];

65

F800GS_L 13



RESULTADOS:

Observando la simulación se puede decir que con las entradas que se usaron la moto está parada y en equilibrio durante los 30 segundos que dura la simulación. Sin embargo la práctica nos indica que esta situación, aunque se puede presentar, no se mantendría durante tanto tiempo.


x0=[Vr,theta1,theta]=[0,0,0];

66

F800GS_L 14

Se simula la respuesta del sistema “F800GSPS”a un vector de entrada u0 determinado. Con base en las señales obtenidas en las simulaciones se determina el vector de estados x0. Para esto se evalúan las señales que hacen parte de dicho vector en un tiempo t = 12 segundos. Wsigma y Fr son entradas paso. El vector de entrada u0 para esta simulación es:

u0=[Fr,Wsigma]=[86,pi/45]

RESULTADOS:

Observando la simulación se puede decir que con las entradas que se usaron la moto pierde el equilibrio poco tiempo después de que Wsigma actúa.

Los vectores de estado y de entrada obtenidos a los 12 segundos de simulación son los siguientes:

x0=[Vr,theta1,theta]=[66.67,0.77,0.35];

67

F800GS_L 15

Se simula la respuesta del sistema “F800GSPS”a un vector de entrada u0 determinado. Con base en las señales obtenidas en las simulaciones se determina el vector de estados x0. Para esto se evalúan las señales que hacen parte de dicho vector en un tiempo t = 12 segundos. Wsigma y Fr son entradas paso. El vector de entrada u0 para esta simulación es :


RESULTADOS:



x0=[Vr,theta1,theta]=[48,0.13,0.065];

68

F800GS_L 16


u0=[Fr,Wsigma]=[86,pi/180];

RESULTADOS:



x0=[Vr,theta1,theta]=[28.1,0.0425,0.028];

69

F800GS_L 17



RESULTADOS:



x0=[Vr,theta1,theta]=[25.24,0.0192,0.0133];

70

F800GS_L 18



RESULTADOS:



x0=[Vr,theta1,theta]=[24.62,0.0094,0.0066];

71

F800GS_L 19


u0=[Fr,Wsigma]=[86,pi/1440];

RESULTADOS:



x0=[Vr,theta1,theta]=[24.48,0.0047,0.0033];

72

F800GS_L 20


u0=[Fr,Wsigma]=[86,pi/2880];

RESULTADOS:



x0=[Vr,theta1,theta]=[24.44,0.0023,0.0016];

73

F800GS_L 21



RESULTADOS:



x0=[Vr,theta1,theta]=[24.42,0.0011,0.0008];

74

ANEXO 3

1. Modelo_Moto

Este es el modelo que se construyó en Simulink™ a partir del modelo que se tomó de la literatura. Ya incluye los saturadores para limitar las variables y es de dos entradas y dos salidas. Las imágenes que vienen a continuación muestran el modelo desde lo más general hasta lo más específico.

Figura A3 - 1

La Figura A2 –2 muestra los bloques que están dentro del bloque Modelo_Moto. En esta figura se puede apreciar claramente que � es la integral de >?, relación que es una parte importante de la definición del modelo. También se pueden ver los saturadores que limitan a #, a -. y a �. Nótese que en este diagrama es claro cómo la entrada >? afecta a todos los bloques de la Figura A2 – 2, ya sea directamente o a través de �.

Figura A3 - 2

75

En la Fgura A3 -3 se muestra una de las ventanas donde se configuración el saturador. Este saturador se usó para limitar a σ.

Figura A3 - 3

La Figura A3 – 4 muestra el bloque NL1 por dentro:

76

Figura A3 - 4

La Figura A3 – 5 muestra el bloque H1 por dentro:

Figura A3 - 5


Figura A3 - 6

77


Figura A3 - 7


Figura A3 - 8

78


Figura A3 - 9

79

2. F800GSPS

A continuación se presenta el modelo simplificado F800GSPS. Con este modelo se llevó a cabo la mayoría de este trabajo. En este modelo ya están incluidos los parámetros que se tomaron de la moto BMW F 800 GS.

Figura A3 - 10

La Figura A3 – 11 muestra el bloque F800GSPS por dentro:

Figura A3 - 11

80


Figura A3 - 12


Figura A3 - 13

81


Figura A3 - 14


Figura A3 - 15

82


Figura A3 - 16


Figura A3 - 17

83

ANEXO 4

1. FS800GS_L1

clear all

%Se definen los valores de los vectores de estado x0 y de entrada u0

%alrededor de los cuales se quiere linealizar el modelo.

Vr = 64.2; theta1 = 0; theta = 0;

Fr = 226; Wsigma = 0;

x0 = [Vr,theta1,theta];

u0 = [Fr,Wsigma];

%% 4.2 Linealización

%Se linealiza el modelo alrededor de los puntos definidos obteniendo como

%resultado su representación en espacio de estados.

[A1,B1,C1,D1] = linmod('F800GSPS',x0,u0);

sys1 = ss(A1,B1,C1,D1);

F_800GS_L1 = pck(A1,B1,C1,D1); %F_800GS_L1 es el modelo linealizado.

%% 4.3.1 Escalización

%Se definen los factores de escalización y se escaliza el modelo en

%espacio de estados obtenido en el paso anterior.

De = mdiag(66.67,pi,pi/2); Du = mdiag(226,pi/45); Dd = mdiag(0,0,pi/9); Dr =

mdiag(66.67,pi,pi/4);

Gs1=mmult(minv(De),F_800GS_L1,Du); %Gs1 es el modelo escalizado.

%% 4.3.2 Realización Mínima

%Se obtiene la realización mínima del modelo escalizado.

[A1L,B1L,C1L,D1L]=unpck(Gs1);

Moto1=ss(A1L,B1L,C1L,D1L);

[mini1,u]=minreal(Moto1);

if isequal(mini1,Moto1)

disp('Gs1 ya es mínima')

end

minGs1 = pck(mini1.a,mini1.b,mini1.c,mini1.d);

%% 4.3.3 Controlabilidad Funcional

%Se genera un vector omega de 3000 puntos espaciado logarítmicamente entre

%0.01 y 1000. Luego se obtiene la respuesta compleja en frecuencia Gf de

%Gs a un vector omega de puntos de frecuencia. Por último se obtiene el

%SVD de la matriz Gf.

omega=logspace(-2,1,300);

Gf1=frsp(minGs1,omega);

[u1,Gfs1,v1]=vsvd(Gf1);

figure;

grid on

axis square

hold on

for i=1:300

Gf1=frsp(minGs1,omega(i));


plot3(u1(1,3),u1(2,3),u1(3,3));

end

figure

hold on

grid on

84

title('Dirección sin ganancia U1(1,3)')

xlabel('rad/s')

ylabel('U1(1,3)')

for i=1:300



plot(omega(i),u1(1,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U1(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U1(3,3)')

for i=1:300




end

%% 4.3.4 Obtención de los polos del sistema

[P, Yp, Xp, Sp] = opde(minGs1);

%% 4.3.5 Obtención de los ceros del sistema

[Z, Yz, Xz, Sz] = ozde(minGs1);

%% 4.3.6 Matriz RGA

Go=frsp(minGs1,0.000001);

RGADC = veval('.*',Go,vpinv(vtp(Go)));

omega1=logspace(-1,3,300);

for i=1:300

minGf1=frsp(minGs1,omega1(i)); %Matriz que representa la

respuesta en frecuencia del sistema escalizado en frecuencia omega(i)

RGA = veval('.*',minGf1,vpinv(vtp(minGf1))); %Saca el RGA en w

RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;

title('RGA(1,1) vs w')

hold on

for i=1:300

plot(omega1(i),RGAM(1,i),'-.r')

end

figure;

85


hold on

for i=1:300

plot(omega1(i),RGAM(2,i),'-.g')

end

figure;


hold on

for i=1:300

plot(omega1(i),RGAM(3,i),'-.b')

end

figure;


hold on

for i=1:300

plot(omega1(i),RGAM(4,i),'-.c')

end

figure;


hold on

for i=1:300

plot(omega1(i),RGAM(5,i),'-.m')

end

figure;


hold on

for i=1:300

plot(omega1(i),RGAM(6,i),'-.k')

end

%% 4.3.7 SVD y respuesta en frecuencia de los valores singulares.

omega1 = logspace(-2,1,300);

Gf171 = frsp(minGs1,omega1);

[u171,Gfs171,v171] = vsvd(Gf171);

figure;

vplot('liv,lm',Gfs171),grid,title('SVD'),xlabel('Rad/s'),ylabel('Ganancia'); %Plot

SVD

2. FS800GS_L2

clear all



Fr = 206; Wsigma = 0;

Vr = 58.5; theta1 = 0; theta = 0;

u0 = [Fr,Wsigma];






sys2 = ss(A2,B2,C2,D2);







86

Gs2 = mmult(minv(De),F_800GS_L2,Du); %Gs2 es el modelo escalizado.



[A2L,B2L,C2L,D2L] = unpck(Gs2);

Moto2 = ss(A2L,B2L,C2L,D2L);

[mini2,u] = minreal(Moto2);



end







omega = logspace(-2,1,300);

Gf2 = frsp(minGs2,omega);

[u2,Gfs2,v2] = vsvd(Gf2);

figure;

grid on

axis square

hold on

for i=1:300



plot3(u2(1,3),u2(2,3),u2(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U2(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U2(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U2(3,3)')

for i=1:300


87



end





%% 4.3.6 Matriz RGA

Go = frsp(minGs2,1);



for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

88




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

3. FS800GS_L3

clear all



Fr = 186; Wsigma = 0;

Vr = 52.8; theta1 = 0; theta = 0;

u0 = [Fr,Wsigma];






sys3 = ss(A3,B3,C3,D3);









%Se obtiene la realización mínima del modelo lineal escalizado.






end










figure;

grid on

axis square

hold on

for i=1:300

89



plot3(u3(1,3),u3(2,3),u3(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U3(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U3(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U3(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA




for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

90

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

4. FS800GS_L4

clear all



Fr = 166; Wsigma = 0;

Vr = 47.15; theta1 = 0; theta = 0;

u0 = [Fr,Wsigma];






91

sys4 = ss(A4,B4,C4,D4);





De = mdiag(66.67,pi,pi/2); Du = mdiag(226,pi/45); Dd = mdiag(0,pi/9,0); Dr =





[A4L,B4L,C4L,D4L]=unpck(Gs4);

Moto4=ss(A4L,B4L,C4L,D4L);

[mini4,u]=minreal(Moto4);



end










figure;

grid on

axis square

hold on

for i=1:300



plot3(u4(1,3),u4(2,3),u4(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U4(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U4(2,3)')

for i=1:300




92

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U4(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA




for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


93

hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

5. FS800GS_L5

clear all



Fr = 146; Wsigma = 0;

Vr = 41.45; theta1 = 0; theta = 0;

u0 = [Fr,Wsigma];






sys5 = ss(A5,B5,C5,D5);















end






94





figure;

grid on

axis square

hold on

for i=1:300



plot3(u5(1,3),u5(2,3),u5(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U5(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U5(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U5(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA




for i=1:300

95




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

6. FS800GS_L6

clear all



Fr = 126; Wsigma = 0;

Vr = 35.78; theta1 = 0; theta = 0;

96

u0 = [Fr,Wsigma];






sys6 = ss(A6,B6,C6,D6);









%Se obtiene la realización mínima del modelo lineal escalizado.






end










figure;

grid on

axis square

hold on

for i=1:300



plot3(u6(1,3),u6(2,3),u6(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U6(1,3)')

for i=1:300




end

figure

hold on

97

grid on


xlabel('rad/s')

ylabel('U6(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U6(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA




for i=1:300

minGf6 = frsp(minGs6,omega1(i)); %Matriz que representa la



RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;

98


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

7. FS800GS_L7

clear all



Fr = 106; Wsigma = 0;

Vr = 30.1; theta1 = 0; theta = 0;

u0 = [Fr,Wsigma];






sys7 = ss(A7,B7,C7,D7);















end

99










figure;

grid on

axis square

hold on

for i=1:300



plot3(u7(1,3),u7(2,3),u7(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U7(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U7(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U7(3,3)')

for i=1:300




end





100

%% 4.3.6 Matriz RGA




for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

8. FS800GS_L8

101

clear all



Fr = 86; Wsigma = 0;

Vr = 24.42; theta1 = 0; theta = 0;

u0 = [Fr,Wsigma];






sys8 = ss(A8,B8,C8,D8);















end










figure;

grid on

axis square

hold on

for i=1:300



plot3(u8(1,3),u8(2,3),u8(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U8(1,3)')

for i=1:300


102



end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U8(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U8(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA




for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;

103


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

9. FS800GS_L9

clear all




Vr = 18.74; theta1 = 0; theta = 0;

u0 = [Fr,Wsigma];






sys9 = ss(A9,B9,C9,D9);











104





end










figure;

grid on

axis square

hold on

for i=1:300



plot3(u9(1,3),u9(2,3),u9(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U9(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U9(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U9(3,3)')

for i=1:300




end


105




%% 4.3.6 Matriz RGA




for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




106

[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

10. FS800GS_L10

clear all




Vr = 13.06; theta1 = 0; theta = 0;

u0 = [Fr,Wsigma];






sys10 = ss(A10,B10,C10,D10);















end









[u10,Gfs10,v10] = vsvd(Gf10);

figure;

grid on

axis square

hold on

for i=1:300



plot3(u10(1,3),u10(2,3),u10(3,3));

end

figure

107

hold on

grid on


xlabel('rad/s')

ylabel('U10(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U10(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U10(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA




for i=1:300

minGf10 = frsp(minGs10,omega1(i)); %Matriz que representa

la respuesta en frecuencia del sistema escalizado en frecuencia omega(i)


RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

108

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

11. FS800GS_L11

clear all




Vr = 7.385; theta1 = 0; theta = 0;

u0 = [Fr,Wsigma];






sys11 = ss(A11,B11,C11,D11);





109











end










figure;

grid on

axis square

hold on

for i=1:300



plot3(u11(1,3),u11(2,3),u11(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U11(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U11(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U11(3,3)')

110

for i=1:300




end





%% 4.3.6 Matriz RGA




for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300

111


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

12. FS800GS_L12

clear all



Fr = 6; Wsigma = 0;

Vr = 1.704; theta1 = 0; theta = 0;

u0 = [Fr,Wsigma];






sys12 = ss(A12,B12,C12,D12);















end









[u12,Gfs12,v12] = vsvd(Gf12);

figure;

grid on

axis square

112

hold on

for i=1:300



plot3(u12(1,3),u12(2,3),u12(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U12(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U12(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U12(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA




for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

113

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

13. FS800GS_L13

14. FS800GS_L14

clear all



Fr = 86; Wsigma = pi/45;

Vr = 66.67; theta1 = 0.77; theta = 0.35;

u0 = [Fr,Wsigma];


114





sys14 = ss(A14,B14,C14,D14);















end









[u14,Gfs14,v14] = vsvd(Gf14);

figure;

grid on

axis square

hold on

for i=1:300

Gf14=frsp(Gs14,omega(i));


plot3(u14(1,3),u14(2,3),u14(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U14(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

115

ylabel('U14(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U14(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA




for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300

116


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

15. FS800GS_L15

clear all




Vr = 48; theta1 = 0.13; theta = 0.65;

u0 = [Fr,Wsigma];






sys15 = ss(A15,B15,C15,D15);















end



117







[u15,Gfs15,v15] = vsvd(Gf15);

figure;

grid on

axis square

hold on

for i=1:300



plot3(u15(1,3),u15(2,3),u15(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U15(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U15(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U15(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA


118



for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

16. FS800GS_L16

clear all

119




Vr = 28.1; theta1 = 0.0425; theta = 0.028;

u0 = [Fr,Wsigma];






sys16 = ss(A16,B16,C16,D16);












[mini16,u] = minreal(Moto16,1.45e-11);



end









[u16,Gfs16,v16] = vsvd(Gf16);

figure;

grid on

axis square

hold on

for i=1:300



plot3(u16(1,3),u16(2,3),u16(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U16(1,3)')

for i=1:300




120

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U16(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U16(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA




for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

121

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

17. FS800GS_L17

clear all




Vr = 25.24; theta1 = 0.0192; theta = 0.0133;

u0 = [Fr,Wsigma];


%% 3.2. Linealización




sys17 = ss(A17,B17,C17,D17);


%% 3.3.1. Escalización






%% 3.3.2. Realización Mínima





122



end


%% 3.3.3. Controlabilidad Funcional



%Gs para un vector omega de puntos de frecuencia. Por último se obtiene el

%SVD de la matriz Gf17.



[u17,Gfs17,v17] = vsvd(Gf17);

figure;

grid on

axis square

hold on

for i=1:300



plot3(u17(1,3),u17(2,3),u17(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U17(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U17(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U17(3,3)')

for i=1:300




end

%% 3.3.4. Obtención de los polos del sistema

123


%% 3.3.5. Obtención de los ceros del sistema


%% Matriz RGA

Go=frsp(minGs17,0,000001);



for i=1:300

minGfs17=frsp(minGs17,omega1(i)); %Matriz que representa

la respuesta en frecuencia del sistema escalizado en frecuencia omega(i)

RGA = veval('.*',minGfs17,vpinv(vtp(minGfs17))); %Saca el RGA en w

RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

%% 4.3. SVD y respuesta en frecuencia de los valores singulares.



124

[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

18. FS800GS_L18

clear all




Vr = 24.62; theta1 = 0.0094; theta = 0.0066;

u0 = [Fr,Wsigma];






sys18 = ss(A18,B18,C18,D18);















end









[u18,Gfs18,v18] = vsvd(Gf18);

figure;

grid on

axis square

hold on

for i=1:300



plot3(u18(1,3),u18(2,3),u18(3,3));

end

figure

125

hold on

grid on


xlabel('rad/s')

ylabel('U18(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U18(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U18(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA

Go = frsp(minGs18,0,000001);



for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

126

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

19. FS800GS_L19

clear all



Fr = 86; Wsigma = pi/1440;

Vr = 66.67; theta1 = 0.77; theta = 0.35;

u0 = [Fr,Wsigma];






sys19 = ss(A19,B19,C19,D19);





127











end









[u19,Gfs19,v19] = vsvd(Gf19);

figure;

grid on

axis square

hold on

for i=1:300



plot3(u19(1,3),u19(2,3),u19(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U19(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U19(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U19(3,3)')

128

for i=1:300




end





%% 4.3.6 Matriz RGA

Go = frsp(minGs19,0.000001);



for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300

129


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

20. FS800GS_L20

clear all



Fr = 86; Wsigma = pi/2880;

Vr = 24.44; theta1 = 0.0023; theta = 0.0016;

u0 = [Fr,Wsigma];






sys20 = ss(A20,B20,C20,D20);















end









[u20,Gfs20,v20] = vsvd(Gf20);

figure;

grid on

axis square

130

hold on

for i=1:300



plot3(u20(1,3),u20(2,3),u20(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U20(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U20(2,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U20(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA

Go = frsp(minGs20,0.000001);



for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

131

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end




[u171,Gfs171,v171] = vsvd(Gf171);

figure;


SVD

21. FS800GS_L21

clear all



Fr = 86; Wsigma = pi/5760;

Vr = 24.42; theta1 = 0.0011; theta = 0.0008;

u0 = [Fr,Wsigma];






132

sys21 = ss(A21,B21,C21,D21);















end









[u21,Gfs21,v21] = vsvd(Gf21);

figure;

grid on

axis square

hold on

for i=1:300



plot3(u21(1,3),u21(2,3),u21(3,3));

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U21(1,3)')

for i=1:300




end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U21(2,3)')

for i=1:300




133

end

figure

hold on

grid on


xlabel('rad/s')

ylabel('U21(3,3)')

for i=1:300




end





%% 4.3.6 Matriz RGA

Go = frsp(minGs21,0.000001);



for i=1:300




RGAM(1,i)=RGA(1,1);

RGAM(2,i)=RGA(2,1);

RGAM(3,i)=RGA(3,1);

RGAM(4,i)=RGA(1,2);

RGAM(5,i)=RGA(2,2);

RGAM(6,i)=RGA(3,2);

end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

for i=1:300


end

figure;


hold on

134

for i=1:300


end

figure;


hold on

for i=1:300


end




[u211,Gfs211,v211] = vsvd(Gf211);

figure;


SVD

135

ANEXO 5

A continuación el código de MATLAB® de las funciones opde y ozde creado por Kjetil Havre para [1], el libro de control multivariable escrito por Sigurd Skogestad e Ian Postlethwaite.

1. Opde

Función para encontrar los polos y direcciones de un sistema representado en espacio de estados.

% [Po, Yp, Xpo, Spo] = opde(G,epp,RF)

%

% The ``Output-Pole-Direction-Eigenvalue'' opde-function computes the

% output pole directions through the use of egenvalue decomposition.

% NOTE: Avoids the change in ordering of eigenvalues that may occur if

% the input and output directions are computed separately.

%

% Outputs: Po - Vector containing the poles.

% Yp - Matrix containing the output pole vectors/directions (normalized =

length 1),

% column i corresponds to pole Po(i).

% Xpo - Output state directions = right eigenvectors of A (non-

normalized).

% NOTE: These are scaled such that pole vectors = pole directions

% Spo - Scalars used to scale the state directions (eigenvectors).

%

% Inputs: G - System matrix on mutools format: G=pck(A,B,C,D)

% epp - (optional) Tolerance, default value is 1e-12.

% It is used to check the norm of the output directions.

% If the norm of output direction is less than EPP,

% zero is stored in the corresponding output direction.

% RF - (optional) If FLG_RF==1 then the first element of

% the output direction is made real.

%

% Written by: Kjetil Havre, 12/9-1996, e-mail: [email protected]

% Requires subroutines: zerm.m and sorte.m

%

% See also: IPDE, PDFST, IZDE, OZDE, ZDSVD and PDSVD.

%

% Reference: Havre K. and S. Skogestad, 1996, ``Effect of RHP Zeros and

% Poles on Performance in Multivariable Systems''.

%

function [P,Y,X,S] = opde(G,epp,rf);

P = []; Y = []; X = [];

[mt,ny,nu, nx] = minfo(G);

S = ones(nx,1);

if strcmp(mt, 'syst')==0

disp( 'System matrix G is required, usage: [] = opde(G,epp)' );

return

end

if nargin < 3

rf = 0;

end

if nargin < 2

epp = 1e-12;

end

[A, B, C, D] = unpck(G);

if norm(D) > 1/epp

disp('Warning: Directions may be inaccurate due to large effect from D.');

end

136

[V, D] = eig(A,'nobalance'); P = diag(D);

[P,Is] = sorte(P); V = V(:,Is);

Hm = A*V-V*diag(P);

if norm(Hm) > 100*epp

disp('Warning: inaccurate eigenvalue computation.')

Hm = Hm;

end

Npd = max(size(P));

for i=1:Npd

X(:,i) = V(:,i);

if abs(P(i)-real(P(i))) < epp

P(i) = real(P(i));

if isreal(X(:,i)) == 0 % Rotate eigenvector.

inz = find(abs(X(:,i))>epp);

if isempty(inz) == 0 % Rotate eigenvector.

vh = angle(X(inz(1),i));

X(:,i)=X(:,i)*exp(-sqrt(-1)*vh);

end

X(:,i) = zerm(X(:,i),epp/10);

end

end

yph = C*X(:,i); nyph = norm(yph);

if nyph > epp

Y(:,i) = yph/nyph;

X(:,i) = X(:,i)/nyph;

S(i,1) = 1/nyph;

if (isreal(Y(:,i))==0) & (rf==1) % Rotate output direction.

inz = find(abs(Y(:,i))>epp);

if isempty(inz) == 0

vh = angle(Y(inz(1),i));

Y(:,i) = Y(:,i)*exp(-sqrt(-1)*vh);

X(:,i)=X(:,i)*exp(-sqrt(-1)*vh);

S(i,1) = S(i,1)*exp(-sqrt(-1)*vh);

end

end

else

Y(:,i) = zeros(ny,1); % Norm is small unobservable.

S(i,1) = 0;

end

end

P = zerm(P,epp/10); X = zerm(X,epp/10); Y = zerm(Y,epp/10);

en

2. Ozde

Función para encontrar los ceros y direcciones de un sistema representado en espacio de estados.

% function [Z, Y, X] = ozde(G,EPP,FLG_RF)

%

% Inputs: G - system matrix in mu-tools format.

% EPP - tolerance, see below, default value 1e-12.

% FLG_RF - If FLG_RF==1 then the first element of

% the output direction is made real.

% Outputs: Z - zeros.

% Y - output zero directions, stored as column vectors.

% X - state directions, stored as column vectors.

% S - S(i) the factor used to normalize output

% direction Y(:,i).

137

% Each column X(:,i) and Y(:,i) corresponds to the zero

% Z(i) and the scaler S(i).

%

% This Output-Zero-Direction-"through genaralized Eigenvalue" decomposition

% OZDE function is a modification of the szeros function contained in mu

% - toolbox. The modification consists of returning the output zero

% directions and the state zero directions in addition to the zeros.

% The output zero directions are defined as:

% y'*G(z) = 0,

% where s = z is a zero of G(s) and ' denotes conjugate transposed.

% This is done by solving the generalized eigenvalue problem of the

% transposed system:

%

% | x' | y' | * | A-Iz | B | = | 0 | 0 |

% |------------|

% | C | D |

%

% are solved by genaralized eigenvalues:

%

% | A"-Iz | C" | * | xi | = | 0 |

% |--------------| |----| |---|

% | B" | D" | | yi | | 0 |

%

% x = conj(xi); y = conj(yi);

%

% The prime ' denotes the conjugate transposed and " denotes the ordinary

% transposed.

%

% OZDE finds the transmission zeros z of a SYSTEM matrix. Occasionally,

% large zeros are included which may actually be at infinity. Solving for

% the transmission zeros of a system involves two generalized eigenvalue

% problems. EPP (optional) defines if the difference between two generalized

% eigenvalues is small. OZDE also finds the output y and the state x

% directions of the zeros.

%

% The output zero directions are stored as column vectors in Y, and the y's

% are normalized so that Y'*Y = I. The state zero directions are stored as

% columns in X. The degree of freedom to normalize the generalized eigen-

% vectors is used to normalize the y part of the vectors. So the length of

% x is not equal to one. Each column in Y and X corresponds to the element

% in Z with same place.

%

% For systems with more outputs than inputs the output zero direction

% is not a complete basis for the left nullspace of G(z).

% Zeros with multiplicity greater than one (rare cases which may

% occure in non-minimal realizations), may (not sure) cause wrong

% directions.

%

% Comments, corrections and malfunctions, can be e-mailed to:

% [email protected] or [email protected]

%

% See also: EIG, SZEROS, IZDE and SPOLES.

% Algorithm based on Laub & Moore 1978 paper, Automatica

%

% Note that when the number of outputs is larger than the number of inputs,

% the output zero direction is not complete. A bit clearer: If z is a zero

% of a non-square plant with number of outputs greater than number of inputs,

% the zero direction is not a line but a surface. As an example, consider

% G(s) with dimensions 3x2 (3 rows and 2 columns). Let s=z be a zero of

% G(s) such that Yz'*G(z) = [0 0]; Since s=z is a zero then the rank of

% G(z) has to be less than the normal rank of G(s), which at maximum can

% be 2. This implies that rank of G(z) must be less than 2.

% y is element in the three dimensional field of real numbers. Since the

% rank of G(z) is maximum one the zero direction is a actually a subspace

138

% in this three dimensional field of real numbers given by two basis

% vectors. Since this function only gives one zero direction for a given

% zero this direction does not describe the zero space on the output

% completely. Two basis vectors are requiered.

%

% Modification for square systems was made by: Kjetil Havre 24/4-1995.

% Modification for non square systems was made by: Kjetil Havre 2/5 -1995.

% Including the state zero dircetions was made by: Kjetil Havre 15/5-1995.

% Modified so that first element of U(:,i) is real: Kjetil Havre 3/2-1996.

% Several modifications made by: Kjetil Havre 5/1-1997.

% - Added optinoal flag RF,

% - Changed the default value of EPP to 1e-12.

% - Function also return the scalers S used to normalize

% the input directions.

% - The zeros ar sorted in ascending order.

%

function [Z, Y, X,S] = ozde(sys,epp,rf)

if nargin < 3

rf = 0;

end

if nargin < 2

epp = 1e-12;

end

if nargin < 1

disp('usage: szeros(sys)')

return

end

[mtype,ny,nu,nx] = minfo(sys);

if mtype == 'syst'

[a,b,c,d] = unpck(sys);

if nx == 0

disp('SYSTEM has no states')

end

sysu = [a b; c d];

sysu = vtp( sysu ); % Transpose system.

% find generalized eigenvalues of a square system matrix

if ny == nu

x = zeros(nx+nu,nx+nu);

x(1:nx,1:nx) = eye(nx);

[vech, ev] = eig(sysu, x);

% Add something here to check for not a number or infinite.

% remove corresponding vectors.

% Also remove top part corresponding to the states

% and normalize bottom part. KH 21/4 - 1995.

z = diag(ev); % Extract the eigenvalues.

kc=0; % Counter for eigenvalues.

for k=1:max(size(z)),

logic = ~isnan(z(k)) & isfinite(z(k));

if logic

kc= kc+1;

Z(kc,1) = z(k);

vech2(:,kc) = vech(:,k);

end

end

Z = zerm(Z,epp/10);

% Split in x and y.

vx = vech2(1:nx, : );

139

vy = vech2(nx+1:nx+ny,:);

% Normalize columns.

[nvr, nvc] = size( vy );

for i=1:nvc,

nrmy = norm(vy(:,i));

if nrmy > 1000*epp

vx(:,i) = vx(:,i)/nrmy;

vy(:,i) = vy(:,i)/nrmy;

S(i,1) = 1/nrmy;

else

vy(:,i) = zeros(ny,1);

S(i,1) = 0;

end

rfi = rf | isreal(Z(i));

if rfi == 1

Inz = find( abs(vy(:,i) ) > 1000*epp );

if isempty(Inz) == 0

X(:,i) = conj(vx(:,i) * exp( -angle(vy(Inz(1),i))*sqrt(-1) ) );

Y(:,i) = conj(vy(:,i) * exp( -angle(vy(Inz(1),i))*sqrt(-1) ) );

S(i,i) = conj(S(i,1) * exp( -angle(vy(Inz(1),i))*sqrt(-1) ) );

else

X(:,i) = conj(vx(:,i));

Y(:,i) = conj(vy(:,i));

end

if isreal(Z(i))

X(:,i) = zerm(X(:,i),epp/10);

Y(:,i) = zerm(Y(:,i),epp/10);

end

else

X(:,i) = conj(vx(:,i));

Y(:,i) = conj(vy(:,i));

end

end

% Find finit zeros.

If = find(abs(Z)< 1/epp );

Z = Z(If);

Y = Y(:,If);

X = X(:,If);

S = S(If);

% Sort the zeros

[Z,Is] = sort(Z);

Y = Y(:,Is);

X = X(:,Is);

S = S(Is);

else % Non-square systems

nrm = norm(sysu,1);

if nu > ny

x1 = [ sysu nrm*(rand(nx+nu,nu-ny)-.5)];

x2 = [ sysu nrm*(rand(nx+nu,nu-ny)-.5)];

else

x1 = [ sysu; nrm*(rand(ny-nu,nx+ny)-.5)];

x2 = [ sysu; nrm*(rand(ny-nu,nx+ny)-.5)];

end

[x]= zeros(size(x1));

x(1:nx,1:nx) = eye(nx);

[v1h z1h] = eig(x1,x); % Compute the genaralized eigenvalues

[v2h z2h] = eig(x2,x); % for the two augumented systems.

z1h2 = diag( z1h );

z2h2 = diag( z2h );

z2 = z2h2(~isnan(z2h2) & isfinite(z2h2));

140

kc=0; % Counter for eigenvalues.

for k=1:max(size(z1h2)),

logic = ~isnan(z1h2(k)) & isfinite(z1h2(k));

if logic

kc= kc+1;

z1(kc,1) = z1h2(k);

vech2(:,kc) = v1h(:,k);

end

end

nz = length(z1);

vech3 = [];

Z = [];

for i=1:nz,

if any(abs(z1(i)-z2) < nrm*sqrt(epp))

Z = [Z; z1(i)];

vech3 = [vech3 vech2(:,i)];

end

end

if isempty(vech3)

Z = []; Y = []; X = [];

return;

end

Z = zerm(Z,epp/10);

% Split columns in x and y.

vx = vech3(1:nx,:);

vy = vech3(nx+1:nx+ny,:);

% Normalize columns.

[nvr, nvc] = size( vy );

for i=1:nvc,

nrmy = norm(vy(:,i));

if nrmy > 1000*epp

vx(:,i) = vx(:,i)/nrmy;

vy(:,i) = vy(:,i)/nrmy;

S(i,1) = 1/nrmy;

else

vy(:,i) = zeros(ny,1);

S(i,1) = 0;

end

rfi = rf | isreal(Z(i));

if rfi == 1

Inz = find( abs(vy(:,i) ) > 1000*epp );

if isempty(Inz) == 0

X(:,i) = conj(vx(:,i) * exp( -angle(vy(Inz(1),i))*sqrt(-1) ) );

Y(:,i) = conj(vy(:,i) * exp( -angle(vy(Inz(1),i))*sqrt(-1) ) );

S(i,1) = conj(S(i,1) * exp( -angle(vy(Inz(1),i))*sqrt(-1) ) );

else

X(:,i) = conj(vx(:,i));

Y(:,i) = conj(vy(:,i));

end

if isreal(Z(i))

X(:,i) = zerm(X(:,i),epp/10);

Y(:,i) = zerm(Y(:,i),epp/10);

end

else

X(:,i) = conj(vx(:,i));

Y(:,i) = conj(vy(:,i));

end

end

% Find finit zeros.

If = find(abs(Z)< 1/epp );

Z = Z(If);

Y = Y(:,If);

X = X(:,If);

141

S = S(If);

% Sort the zeros

[Z,Is] = sort(Z);

Y = Y(:,Is);

X = X(:,Is);

S = S(Is);

end

else

error('matrix is not a SYSTEM matrix')

return

end

%

% Copyright MUSYN INC 1991, All Rights Reserved

% Copyright MUSYN INC 1995, All Rights Reserved

%

ANÁLISIS DE CONTROLABILIDAD PARA UN MODELO DINÁMICO DE UNA …

Documents

Transcript of ANÁLISIS DE CONTROLABILIDAD PARA UN MODELO DINÁMICO DE UNA …