UNIVERSIDAD PERUANA DE CIENCIAS APLICADAS

FACULTAD DE INGENIERÍA

ESCUELA DE INGENIERÍA MECATRÓNICA

INGENIERÍA DE CONTROL 3

LABORATORIO CALIFICADO 2

INTEGRANTES

Aquiño Aquino, Herber Giovanni U20201E912

Sanca Condori, Sandro U20211C259
Alcca Ramirez, Luis Maximiliano U20211B184
Altamirano Casas, Ernesto Alonso U201614377

PROFESOR
Carlos Hernan Inga Espinoza

Sección – MS83

14 de Noviembre del 2024

Lima - Perú
 ÍNDICE

INTRODUCCIÓN ........................................................................................................ 4
MODELAMIENTO MATEMÁTICO DEL MOTOR DC EN ESPACIO DE ESTADOS .... 6
1. Ecuaciones Dinámicas......................................................................................... 6
2. Espacio de Estados ............................................................................................. 7
a) Velocidad 𝝎𝒕 .................................................................................................. 7
b) Posición 𝜽𝒕 ..................................................................................................... 7
TIEMPO DE MUESTREO Y MODELO MATEMÁTICO EN TIEMPO DISCRETO....... 8
1. Tiempo de Muestreo ........................................................................................... 8
2. Modelo Discreto ............................................................................................... 10
REQUERIMIENTOS DE DISEÑO DEL CONTROLADOR ......................................... 11
DISEÑO Y SIMULACIÓN DE LOS CONTROLADORES............................................ 11
1. Controlador por Reubicación de Polos ............................................................... 11
a) Velocidad ..................................................................................................... 11
b) Posición ........................................................................................................ 12
2. Controlador por Reubicación de Polos con Acción integrativa ............................. 13
a) Velocidad ..................................................................................................... 13
b) Posición ........................................................................................................ 14
3. Controlador Óptimo ......................................................................................... 17
a) Velocidad ..................................................................................................... 17
b) Posición ........................................................................................................ 19
4. Controlador Óptimo con Acción integrativa ....................................................... 20
a) Velocidad ..................................................................................................... 20
b) Posición ........................................................................................................ 21
INSTRUMENTACIÓN Y DISEÑO DEL OBSERVADOR (FILTRO DE KALMAN) ...... 23
1. Instrumentación Mínima .................................................................................. 23
2. Filtro de Kalman .............................................................................................. 24
FUSIÓN DE SENSORES ............................................................................................ 24
INTEGRACIÓN CONTROLADOR-OBSERVADOR ................................................... 26
DIAGRAMA DE SIMULACIÓN EN TIEMPO DISCRETO ......................................... 27
CONCLUSIONES ...................................................................................................... 27
ANEXOS ................................................................................................................... 28
INTRODUCCIÓN
En este laboratorio, se abordará el diseño y simulación de un sistema de control para un motor de
corriente continua (DC). Este tipo de motor es ampliamente utilizado en aplicaciones industriales
y de automatización debido a su simplicidad y facilidad de control. No obstante, su regulación
precisa requiere el desarrollo de estrategias de control avanzadas que permitan estabilizar su
velocidad y posición ante perturbaciones externas.

El objetivo principal de este laboratorio es diseñar y evaluar diferentes controladores digitales

mediante retroalimentación de estados, asegurando que el sistema responda de manera eficiente y
estable. Para lograr esto, se va a desarrollar y simular varios métodos de control, tales como el
controlador por reubicación de polos y el controlador óptimo, en sus versiones con y sin acción
integrativa. Adicionalmente, se realizará la implementación de un observador basado en el filtro
de Kalman, que permita estimar el estado del sistema en presencia de ruido, y así evaluar el
desempeño de la fusión de sensores como técnica para mejorar la precisión de las mediciones.

Asimismo, se abordarán los pasos fundamentales del diseño, desde el modelamiento matemático
del motor DC en espacio de estados hasta la implementación del sistema completo en tiempo
discreto. Se seleccionará un tiempo de muestreo adecuado, y se discretizarán las ecuaciones del
modelo para construir un controlador digital robusto. Finalmente, se integrarán el controlador y el
observador en una configuración de lazo cerrado, y se analizarán los resultados obtenidos para
evaluar el desempeño global del sistema.

Los datos que se usará en este laboratorio se presentan en la figura 1, siendo en específico para
nuestro grupo la columna 2 (222049):
Fig. 1: Parámetros motor DC.
MODELAMIENTO MATEMÁTICO DEL MOTOR DC EN ESPACIO DE ESTADOS

1. Ecuaciones Dinámicas
Para modelar un motor de corriente continua (DC), consideramos las ecuaciones que describen la
relación entre el voltaje aplicado, la corriente del motor, el torque, y la velocidad angular. Estas
ecuaciones se basan en dos aspectos fundamentales: la ecuación eléctrica y la ecuación mecánica.

Ecuación Eléctrica

La ecuación que describe el circuito eléctrico del motor DC se basa en la Ley de Kirchhoff para la
tensión:
𝒅
𝑽(𝒕) = 𝑹 ∙ 𝒊(𝒕) + 𝑳 ∙ 𝒊(𝒕) + 𝑬𝒃 (𝒕)
𝒅𝒕
Donde:

𝑽(𝒕): 𝑇𝑒𝑛𝑠𝑖ó𝑛 𝑎𝑝𝑙𝑖𝑐𝑎𝑑𝑎 𝑎𝑙 𝑚𝑜𝑡𝑜𝑟 (𝑉).

𝑹: 𝑅𝑒𝑠𝑖𝑠𝑡𝑒𝑛𝑐𝑖𝑎 𝑑𝑒𝑙 𝑚𝑜𝑡𝑜𝑟 (𝛺).

𝑳: 𝐼𝑛𝑑𝑢𝑐𝑡𝑎𝑛𝑐𝑖𝑎 𝑑𝑒𝑙 𝑚𝑜𝑡𝑜𝑟 (𝐻).

𝒊(𝒕): 𝐶𝑜𝑟𝑟𝑖𝑒𝑛𝑡𝑒 𝑎 𝑡𝑟𝑎𝑣é𝑠 𝑑𝑒𝑙 𝑚𝑜𝑡𝑜𝑟 (𝐴).

𝑬𝒃 (𝒕): Fuerza contraelectromotriz (V), 𝑞𝑢𝑒 𝑒𝑠 𝑝𝑟𝑜𝑝𝑜𝑟𝑐𝑖𝑜𝑛𝑎𝑙 𝑎 𝑙𝑎 𝑣𝑒𝑙𝑜𝑐𝑖𝑑𝑎𝑑 𝑎𝑛𝑔𝑢𝑙𝑎𝑟 𝜔(𝑡).

La fuerza contraelectromotriz 𝐸𝑏 (𝑡) está relacionada con la velocidad angular del motor
𝜔(𝑡) mediante la constante de fuerza contraelectromotriz 𝐾𝑏 :

𝑬𝒃 (𝒕) = 𝑲𝒃 ∙ 𝝎(𝒕)

Ecuación Mecánica

La ecuación de movimiento para la parte mecánica del motor se describe con la Segunda Ley de
Newton para la rotación:

𝑱 ∙ 𝝎̇(𝒕) = 𝑻(𝒕) − 𝑩𝒎 ∙ 𝝎(𝒕)

Donde:

𝑱: 𝑀𝑜𝑚𝑒𝑛𝑡𝑜 𝑑𝑒 𝑖𝑛𝑒𝑟𝑐𝑖𝑎 𝑑𝑒𝑙 𝑚𝑜𝑡𝑜𝑟 (𝑘𝑔 ∙ 𝑚2 ).

𝝎(𝒕): 𝑉𝑒𝑙𝑜𝑐𝑖𝑑𝑎𝑑 𝑎𝑛𝑔𝑢𝑙𝑎𝑟 (𝑟𝑎𝑑/𝑠).

𝑻(𝒕): 𝑇𝑜𝑟𝑞𝑢𝑒 𝑔𝑒𝑛𝑒𝑟𝑎𝑑𝑜 𝑝𝑜𝑟 𝑒𝑙 𝑚𝑜𝑡𝑜𝑟(𝑁 ∙ 𝑚).

𝑩𝒎 : 𝐶𝑜𝑒𝑓𝑖𝑐𝑖𝑒𝑛𝑡𝑒 𝑑𝑒 𝑓𝑟𝑖𝑐𝑐𝑖ó𝑛 𝑣𝑖𝑠𝑐𝑜𝑠𝑎 (𝑁 ∙ 𝑚 ∙ 𝑠/𝑟𝑎𝑑).

El torque 𝑇(𝑡) es proporcional a la corriente 𝑖(𝑡) a través de la constante de torque 𝐾𝑡 :

𝑻(𝒕) = 𝑲𝒕 ∙ 𝒊(𝒕)

También, se sabe que:

𝝎(𝒕) = 𝜽̇(𝒕)

2. Espacio de Estados

a) Velocidad 𝝎(𝒕)
❖ Variables de Estado
𝜔(𝑡): 𝑉𝑒𝑙𝑜𝑐𝑖𝑑𝑎𝑑 𝐴𝑛𝑔𝑢𝑙𝑎𝑟
𝑖(𝑡): 𝐶𝑜𝑟𝑟𝑖𝑒𝑛𝑡𝑒

𝒙𝟏 𝝎(𝒕)
𝒙 = [𝒙 ] = [ ]
𝟐 𝒊(𝒕)

❖ Ecuaciones
𝑲𝒕 𝑩𝒎
𝒙̇ 𝟏 (𝒕) = ∙ 𝒙𝟐 (𝒕) − ∙ 𝒙𝟏 (𝒕)
𝑱 𝑱
𝑹 𝑲𝒃 𝟏
𝒙̇ 𝟐 (𝒕) = − ∙ 𝒙𝟐 (𝒕) − ∙ 𝒙𝟏 (𝒕) + ∙ 𝑽(𝒕)
𝑳 𝑳 𝑳

❖ Representación
𝒙̇ = 𝑨𝒙 + 𝑩𝒖

𝑩𝒎 𝑲𝒕 𝟎
−
𝒙̇ 𝟏 𝑱 𝑱 𝒙𝟏
[ ]= [𝒙 ] + [𝑽(𝒕)]
𝒙̇ 𝟐 𝟐 𝟏
𝑲𝒃 𝑹
−
[ 𝑳 − ] [𝑳 ]
𝑳

b) Posición 𝜽(𝒕)
❖ Variables de Estado
 𝜃(𝑡): 𝑃𝑜𝑠𝑖𝑐𝑖ó𝑛 𝐴𝑛𝑔𝑢𝑙𝑎𝑟
𝜃̇(𝑡): 𝑉𝑒𝑙𝑜𝑐𝑖𝑑𝑎𝑑 𝐴𝑛𝑔𝑢𝑙𝑎𝑟
𝑖(𝑡): 𝐶𝑜𝑟𝑟𝑖𝑒𝑛𝑡𝑒

𝒙𝟏 𝜽(𝒕)
𝒙 = [ 𝟐 ] = [𝜽̇(𝒕)]
𝒙
𝒙𝟑 𝒊(𝒕)

❖ Ecuaciones
𝑲𝒕 𝑩𝒎
𝜽̈(𝒕) = ∙ 𝒊(𝒕) − ∙ 𝜽̇(𝒕)
𝑱 𝑱
𝒅 𝑹 𝑲𝒃 𝟏
𝒊(𝒕) = − ∙ 𝒊(𝒕) − ∙ 𝜽̇(𝒕) + ∙ 𝑽(𝒕)
𝒅𝒕 𝑳 𝑳 𝑳

❖ Representación
𝒙̇ = 𝑨𝒙 + 𝑩𝒖

𝟎 𝟏 𝟎 𝟎

𝒙̇ 𝟏 𝑩𝒎 𝑲𝒕 𝒙𝟏
[𝒙̇ 𝟐 ] = 𝟎 −
𝑱
𝒙
𝑱 [ 𝟐] +
𝟎 [𝑽(𝒕)]
𝒙̇ 𝟑 𝒙𝟑
𝑲𝒃 𝑹 𝟏
[𝟎 −
𝑳
− ]
𝑳 [ 𝑳]

TIEMPO DE MUESTREO Y MODELO MATEMÁTICO EN TIEMPO DISCRETO

1. Tiempo de Muestreo
A partir de la figura 1, los parámetros del motor DC serán los del modelo 222049:

• 𝑅 = 3,18 Ω
• 𝐿 = 0,154*10-3 𝐻
• 𝐽 = 4.22 ∙ 10−7 𝑘𝑔 ∙ 𝑚2
• 𝐾𝑡 = 0.0151 N m / A
• Kb = 0.0151 v/rad*s
• B = 0.0221*10-3 𝑁 ∙ 𝑚 ∙ 𝑠

𝜔(𝑠) 2.3235 ∗ 108

= 2
𝑉(𝑠) 𝑆 + 2.0702 ∗ 104 𝑆 + 4.5899 ∗ 106
Tenemos:

• 𝜔𝑛 = 2142.40519

• 𝜉 = 4.83

• 𝑘 = 50.62

Entonces tenemos que el tiempo de establecimiento es:

4
𝑇𝑠𝑡 = = 3.87 ∗ 10−4 𝑠
ξ ∗ Wn

A partir del tiempo de establecimiento, podemos calcular el tiempo de muestreo:

𝑇𝑠𝑡
𝑇= , 𝑁𝑟 ∈ [25,75] → 𝑇 = (0.55 𝑉 5.16)𝜇𝑠
𝑁𝑟
Como el rango del tiempo de muestre resultó ser bastante pequeño, escogeremos:

𝑇 = 5 𝜇𝑠

𝜃(𝑠) 2.324 ∗ 108

= 2
𝑉(𝑠) 𝑠 + 2.07 ∗ 104 𝑠 + 4.59 ∗ 106

Tenemos:

• 𝜔𝑛 = 2142.428529

• 𝜉 = 4.83

• 𝑘 = 50.63180828

Entonces tenemos que el tiempo de establecimiento es:

4
𝑇𝑠𝑡 = = 3.87 ∗ 10−4 𝑠
ξ ∗ Wn

A partir del tiempo de establecimiento, podemos calcular el tiempo de muestreo:

𝑇𝑠𝑡
𝑇= , 𝑁𝑟 ∈ [25,75] → 𝑇 = (0.155 𝑉 5.16)𝜇𝑠
𝑁𝑟
Como el rango del tiempo de muestreo resultó ser bastante pequeño, escogeremos:

𝑇 = 5 𝜇𝑠

2. Modelo Discreto
Velocidad:

𝑥̇ (𝑡) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡)

𝑦(𝑡) = 𝐶𝑥(𝑡) + 𝐷𝑢(𝑡)

−𝐵 𝐾𝑡
𝐽 𝐽
𝐴=
−𝐾𝑏 −𝑅
[ 𝐿 𝐿 ]
0
𝐵 = [ 1]
𝐿

𝐶 = [ 1 0]

𝑋1 (𝑘 + 1) 0.9997 0.1700 𝑋1 (𝑘) 0.0028

[ ]= [ ][ ]+ [ ] 𝑢(𝑘)
𝑋2 (𝑘 + 1) −0.0005 0.9019 𝑋2 (𝑘) 0.0308

𝑋1 (𝑘)
𝑦(𝑘) = [1 0] [ ]
𝑋2 (𝑘)

Posición:
𝑥̇ (𝑡) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡)

𝑦(𝑡) = 𝐶𝑥(𝑡) + 𝐷𝑢(𝑡)

−𝑅 −𝐾𝑏
0
𝐿 𝐿
𝐴 = 𝐾𝑡 −𝐵
0
𝐽 𝐽
[ 0 1 0]
 1
𝐵 = [𝐿 ]
0
0

𝐶 = [0 1 0]

𝑋1 (𝑘 + 1) 0.9019 −0.0003 0 𝑋1 (𝑘) 0.0308

[𝑋2 (𝑘 + 1)] = [0.1139 0.4296 0] [𝑋2 (𝑘)] + [0.0022] 𝑢(𝑘)
𝑋3 (𝑘 + 1) 0 0 1 𝑋3 (𝑘) 0

𝑋1 (𝑘)
𝑦(𝑘) = [0 1 0] [𝑋2 (𝑘)]
𝑋3 (𝑘)

REQUERIMIENTOS DE DISEÑO DEL CONTROLADOR

Para los requerimientos del controlador tendremos los siguientes:

Tiempo de establecimiento del sistema (𝒕𝒔𝒕: tiempo en el que el sistema alcanza un punto estable
sin fluctuaciones): 2 segundos

Máximo sobre impulso (𝑴𝒑: máximo desvío del primer pico de impulso): 10%

DISEÑO Y SIMULACIÓN DE LOS CONTROLADORES

1. Controlador por Reubicación de Polos

a) Velocidad
PASO 1: Controlabilidad (Co)

𝐶𝑂 = [𝐵 𝐴𝐵 ]

𝑅𝑎𝑛𝑔𝑜(𝐶𝑜 ) = 2

PASO 2: Polinomio Característico Original

|𝑠𝐼 − 𝐴| = 𝑠 2 + 𝑎1 𝑠 + 𝑎2
 PASO 3: Matriz W
𝑎1 1
𝑊=[ ]
1 0

PASO 4: Matriz de transformación (T)

𝑇 = 𝐶𝑜 × 𝑊

PASO 5: Polinomio deseado

𝑠 2 + 𝑚1 𝑠 + 𝑚2 = 0

PASO 6: Matriz de Ganancia (𝑲𝒁 )

𝐾𝑧 = [𝑚2 − 𝑎2 𝑚1 − 𝑎1 ]

PASO 7: Matriz de Ganancia Controlador (K)

𝐾 = 𝐾𝑍 𝑇 −1

PASO 8: Ley de Control

𝑢 = 𝐾(𝑟 − 𝑥)

b) Posición
PASO 1: Controlabilidad (Co)

𝐶𝑂 = [𝐵 𝐴𝐵 𝐴2 𝐵 ]

𝑅𝑎𝑛𝑔𝑜(𝐶𝑜 ) = 3

PASO 2: Polinomio Característico Original

|𝑠𝐼 − 𝐴| = 𝑠 3 + 𝑎1 𝑠 2 + 𝑎2 𝑠 + 𝑎3

PASO 3: Matriz W
 𝑎2 𝑎1 1
𝑊 = [𝑎1 1 0]
1 0 0

PASO 4: Matriz de transformación (T)

𝑇 = 𝐶𝑜 × 𝑊

PASO 5: Polinomio deseado

𝑠 3 + 𝑚1 𝑠 2 + 𝑚2 𝑠 + 𝑚3 = 0

PASO 6: Matriz de Ganancia (𝑲𝒁 )

𝐾𝑧 = [𝑚3 − 𝑎3 𝑚2 − 𝑎 2 𝑚1 − 𝑎1 ]

PASO 7: Matriz de Ganancia Controlador (K)

𝐾 = 𝐾𝑍 𝑇 −1

PASO 8: Ley de Control

𝑢 = 𝐾(𝑟 − 𝑥)

2. Controlador por Reubicación de Polos con Acción integrativa

a) Velocidad
PASO 1: Controlabilidad (Co)

𝐶𝑂 = [𝐵𝑖 𝐴𝑖 𝐵𝑖 𝐴2𝑖 𝐵𝑖 ]

PASO 2: Polinomio Característico Original

|𝑠𝐼 − 𝐴𝑖 | = 𝑠 3 + 𝑎1 𝑠 2 + 𝑎2 𝑠 + 𝑎3

PASO 3: Matriz W
 𝑎2 𝑎1 1
𝑊 = [𝑎1 1 0]
1 0 0

PASO 4: Matriz de transformación (T)

𝑇 = 𝐶𝑜 × 𝑊

PASO 5: Polinomio deseado

𝑠 3 + 𝑚1 𝑠 2 + 𝑚2 𝑠 + 𝑚3 = 0

PASO 6: Matriz de Ganancia (𝑲𝒁 )

𝐾𝑧 = [𝑚3 − 𝑎3 𝑚2 − 𝑎 2 𝑚1 − 𝑎1 ]

PASO 7: Matriz de Ganancia Controlador (K)

𝐾 = 𝐾𝑍 𝑇 −1

PASO 8: Ley de Control

𝒖 = −𝑲(𝟏:𝟐) ∙ 𝒙 + 𝑲(𝟑) ∫ 𝒆 𝒅𝒕

b) Posición
PASO 1: Controlabilidad (Co)

𝐶𝑂 = [𝐵𝑖 𝐴𝑖 𝐵𝑖 𝐴2𝑖 𝐵𝑖 𝐴3𝑖 𝐵𝑖 ]

PASO 2: Polinomio Característico Original

|𝑠𝐼 − 𝐴𝑖 | = 𝑠 4 + 𝑎1 𝑠 3 + 𝑎2 𝑠 2 + 𝑎3 𝑠 + 𝑎4

PASO 3: Matriz W
𝑎3 𝑎2 𝑎1 1
𝑎 𝑎1 1 0
𝑊=[ 2 ]
𝑎1 1 0 0
1 0 0 0
PASO 4: Matriz de transformación (T)

𝑇 = 𝐶𝑜 × 𝑊

PASO 5: Polinomio deseado

𝑠 4 + 𝑚1 𝑠 3 + 𝑚2 𝑠 2 + 𝑚3 𝑠 + 𝑚4 = 0

PASO 6: Matriz de Ganancia (𝑲𝒁 )

𝐾𝑧 = [𝑚4 − 𝑎4 𝑚3 − 𝑎3 𝑚2 − 𝑎 2 𝑚1 − 𝑎1 ]

PASO 7: Matriz de Ganancia Controlador (K)

𝐾 = 𝐾𝑍 𝑇 −1

PASO 8: Ley de Control

𝒖 = −𝑲(𝟏:𝟑) ∙ 𝒙 + 𝑲(𝟒) ∫ 𝒆 𝒅𝒕
 Simulación de controlador por reubicación de polos con y sin AI de velocidad

Fig. 2: Simulación de controlador por reubicación de polos de velocidad.

El gráfico de la simulación nos muestra la diferencia entre el controlador por reubicación de polos
con acción integrativa, y el controlador por reubicación de polos sin acción integrativa, ambos de
la velocidad.
 Simulación de controlador por reubicación de polos con y sin AI de posición

Fig. 3: Simulación de controlador por reubicación de polos de posición.

El gráfico de la simulación nos muestra la diferencia entre el controlador por reubicación de polos
con acción integrativa, y el controlador por reubicación de polos sin acción integrativa, ambos de
la posición. Notamos que la acción integrativa tiene menos sobre impulso que sin acción
integrativa pero los dos se establecen en cierto tiempo.

3. Controlador Óptimo
El diseño de este controlador implica la definición de una función de costo y la solución de la
ecuación de Riccati asociada para obtener la ganancia de retroalimentación del estado K.

a) Velocidad
Metodología de control óptimo
𝜔(𝑡) 𝑥1
𝑥=[ ] = [𝑥 ] , 𝑢 = [𝑣(𝑡)] = [𝑢1 ]
𝑖(𝑡) 2

Tenemos nuestra función de costo:

𝑡𝑓
1 𝑇 1 𝑇
𝐽(𝑥, 𝑢) = 𝑥(𝑡𝑓) 𝑄𝑥(𝑡𝑓) + ∫ (𝑥 𝑄𝑥 + 𝑢𝑇 𝑅𝑢)𝑑𝑡
2 0 2
Donde:

La matriz Q es una matriz de pesos que penaliza los errores en los estados del
sistema.

- Valores grandes de Q: minimiza el error.

- Ajustando valores de Q se puede dar más o menos importancia a los
estados del sistema.

La matriz R es una matriz de pesos que penaliza el esfuerzo de control.

- Valores de R grandes: movimientos más suaves o menos energía

consumida por el actuador.
- Ajustar R para que el controlador no exceda los límites del actuador
para el control.

Para que 𝐽(𝑥, 𝑢) sea convexa y positiva:

𝑄>0

𝑅>0

Tenemos:
𝑞 0
𝑄=[ 1 ]
0 𝑞2

𝑅 = [𝑟1 ]

Entonces:

𝑞1 𝑦 𝑞2 > 0

𝑟1 > 0

Ecuación de Riccati

𝐴𝑇 𝑃 + 𝑃𝐴 − 𝑃𝐵𝑅 −1 𝐵 𝑇 𝑃 + 𝑄 = 0

𝑃 𝑒𝑠 𝑙𝑎 𝑠𝑜𝑙𝑢𝑐𝑖𝑜𝑛 𝑑𝑒 𝑙𝑎 𝑒𝑐𝑢𝑎𝑐𝑖ó𝑛 𝑑𝑒 𝑅𝑖𝑐𝑐𝑎𝑡𝑖.

Y se obtiene en el software Matlab de la siguiente manera:

𝑃 = 𝑎𝑟𝑒(𝐴, 𝐵𝑅 −1 𝐵 𝑇 , 𝑄)
 Ganancia del controlador K:

𝐾 = 𝑅 −1 𝐵𝑇 𝑃

Ley de control:

𝑢 = 𝐾 × (𝑟𝑒𝑓 − 𝑥𝑜𝑝 )

b) Posición
Metodología de control óptimo
𝜃(𝑡) 𝑥1
𝑥 = [𝜃̇ (𝑡)] = [𝑥2 ] , 𝑢 = [𝑣(𝑡)] = [𝑢1 ]
𝑖(𝑡) 𝑥3

Tenemos nuestra función de costo:

𝑡𝑓
1 𝑇 1 𝑇
𝐽(𝑥, 𝑢) = 𝑥(𝑡𝑓) 𝑄𝑥(𝑡𝑓) + ∫ (𝑥 𝑄𝑥 + 𝑢𝑇 𝑅𝑢)𝑑𝑡
2 0 2

Para que 𝐽(𝑥, 𝑢) sea convexa y positiva:

𝑄>0

𝑅>0

Tenemos:
𝑞1 0 0
𝑄 = [0 𝑞2 0]
0 0 𝑞3

𝑅 = [𝑟1 ]

Entonces:

𝑞1 , 𝑞2 𝑦 𝑞3 > 0

𝑟1 > 0

Ecuación de Riccati
 𝐴𝑇 𝑃 + 𝑃𝐴 − 𝑃𝐵𝑅 −1 𝐵 𝑇 𝑃 + 𝑄 = 0

𝑃 𝑒𝑠 𝑙𝑎 𝑠𝑜𝑙𝑢𝑐𝑖𝑜𝑛 𝑑𝑒 𝑙𝑎 𝑒𝑐𝑢𝑎𝑐𝑖ó𝑛 𝑑𝑒 𝑅𝑖𝑐𝑐𝑎𝑡𝑖.

Y se obtiene en el software Matlab de la siguiente manera:

𝑃 = 𝑎𝑟𝑒(𝐴, 𝐵𝑅 −1 𝐵 𝑇 , 𝑄)

Ganancia del controlador K:

𝐾 = 𝑅 −1 𝐵𝑇 𝑃

Ley de control:

𝑢 = 𝐾 × (𝑟𝑒𝑓 − 𝑥𝑜𝑝 )

4. Controlador Óptimo con Acción integrativa

Para el diseño del controlador óptimo con acción integrativa, al igual que con el controlador
por reubicación de polos con acción integrativa se amplían las matrices para su diseño y
obtener las ganancias correspondientes para la ley de control.

a) Velocidad
Metodología de control óptimo con acción integrativa

Matriz ampliada
𝟎 𝑩
𝑨𝒊 = [ 𝑨 𝟎] , 𝑩𝒊 = [ ]
𝟎
𝟏 𝟎 𝟎

Ecuación de Riccati

𝐴𝑇𝑖 𝑃 + 𝑃𝐴𝑖 − 𝑃𝐵𝑖 𝑅 −1𝐵𝑖𝑇 𝑃 + 𝑄𝑖 = 0

𝑃 𝑒𝑠 𝑙𝑎 𝑠𝑜𝑙𝑢𝑐𝑖𝑜𝑛 𝑑𝑒 𝑙𝑎 𝑒𝑐𝑢𝑎𝑐𝑖ó𝑛 𝑑𝑒 𝑅𝑖𝑐𝑐𝑎𝑡𝑖.

Y se obtiene en el software Matlab de la siguiente manera:

𝑃 = 𝑎𝑟𝑒(𝐴𝑖 , 𝐵𝑖 𝑅 −1 𝐵𝑖𝑇 , 𝑄𝑖 )
 Ganancia del controlador Ki:

𝐾𝑖 = 𝑅 −1 𝐵𝑖𝑇 𝑃

Ley de control:

𝑢 = −𝐾𝑖(1:2) 𝑥𝑜𝑝𝑖 + 𝐾𝑖(3) ∫ 𝑒 𝑑𝑡

b) Posición
Metodología de control óptimo con acción integrativa

Matriz ampliada
𝟎
𝑩
𝑨𝒊 = [ 𝑨 𝟎] , 𝑩𝒊 = [ ]
𝟎 𝟎
𝟏𝟎𝟎 𝟎

Ecuación de Riccati

𝐴𝑇𝑖 𝑃 + 𝑃𝐴𝑖 − 𝑃𝐵𝑖 𝑅 −1𝐵𝑖𝑇 𝑃 + 𝑄𝑖 = 0

𝑃 𝑒𝑠 𝑙𝑎 𝑠𝑜𝑙𝑢𝑐𝑖𝑜𝑛 𝑑𝑒 𝑙𝑎 𝑒𝑐𝑢𝑎𝑐𝑖ó𝑛 𝑑𝑒 𝑅𝑖𝑐𝑐𝑎𝑡𝑖.

Y se obtiene en el software Matlab de la siguiente manera:

𝑃 = 𝑎𝑟𝑒(𝐴𝑖 , 𝐵𝑖 𝑅 −1 𝐵𝑖𝑇 , 𝑄𝑖 )

Ganancia del controlador Ki:

𝐾𝑖 = 𝑅 −1 𝐵𝑖𝑇 𝑃

Ley de control:

𝑢 = −𝐾𝑖(1:3) 𝑥𝑜𝑝𝑖 + 𝐾𝑖(4) ∫ 𝑒 𝑑𝑡

 Simulación de controlador optimo con y sin AI de velocidad

Fig. 4: Simulación de controlador optimo de velocidad.

El gráfico de la simulación nos muestra la diferencia entre el controlador optimo con acción
integrativa, y el controlador optimo sin acción integrativa, ambos de la velocidad. Notamos que
sin acción integrativa tiene una respuesta más rápida que con acción integrativa.
 Simulación de controlador optimo con y sin AI de posición

Fig. 5: Simulación de controlador optimo de posición.

El gráfico de la simulación nos muestra la diferencia entre el controlador optimo con acción
integrativa, y el controlador optimo sin acción integrativa, ambos de la posición. Notamos que la
acción integrativa tiene menos tiempo de respuesta que sin acción integrativa.

INSTRUMENTACIÓN Y DISEÑO DEL OBSERVADOR (FILTRO DE KALMAN)

1. Instrumentación Mínima
Sensor de posición: Un encoder rotativo es esencial para medir la posición angular del motor.

Sensor de velocidad: Puede ser derivado del encoder rotativo o utilizar un sensor específico como
un tacómetro si se requiere mayor precisión.
Convertidor A/D y D/A: Para digitalizar las señales analógicas provenientes de los sensores y
para enviar señales analógicas al motor, si es necesario.

Controlador (Microcontrolador o DSP): Un dispositivo capaz de ejecutar el controlador óptimo

y el filtro de Kalman. Un microcontrolador como un STM32 o un microprocesador como un DSP
sería adecuado.

H-Bridge o controlador de motor: Para manejar la dirección y velocidad del motor DC.

Fuente de alimentación: Una adecuada para el motor y los componentes electrónicos.

2. Filtro de Kalman

a) Predicción del estado:

b) Predicción de la covarianza:

c) Ganancia de Kalman:

d) Actualización del estado:

e) Actualización de la covarianza:

FUSIÓN DE SENSORES
Realizamos la fusión de sensores por el método de media ponderada, con un peso asignado a cada
uno. Los sensores que utilizaremos son el encoder, por su precio económico, pero presentando
ruido en bajas velocidades, y el tacómetro, este es menos preciso, pero con una respuesta más
estable a bajas velocidades.

Por la media ponderada, asignamos pesos llamados 𝜔1 y 𝜔2 para las mediciones de cada sensor,
de acuerdo con su precisión relativa. En caso del encoder, consideramos una precisión del 90% y
del tacómetro un 70%, entonces los pesos relativos se hallarían así:
0.9
𝜔1 = = 0.5625
0.9 + 0.7
0.7
𝜔2 = = 0.4375
0.9 + 0.7
Consideramos el encoder como 𝑋1 y al tacómetro como 𝑋2, la estimación combinada es 𝑋̇:

𝑋̇ = 𝜔1 𝑋1 + 𝜔2 𝑋2

Esto pondera la medición del sensor con mayor precisión, reduciendo la influencia de mediciones
menos confiables.

Implementando la fusión de los sensores, obtenemos la siguiente respuesta:

Fig. 6: Simulación de fusión de sensores.

INTEGRACIÓN CONTROLADOR-OBSERVADOR

Fig. 7: Simulación de planta controlada usando la señal original y la señal observada.

Fig. 8: Simulación de planta controlada usando la señal observada.

DIAGRAMA DE SIMULACIÓN EN TIEMPO DISCRETO

Fig. 9: Planta controlada usando la señal de velocidad original y del controlador.

CONCLUSIONES
La implementación de controladores por reubicación de polos y controladores óptimos, con y sin
acción integrativa, permitió mejorar la precisión y estabilidad del sistema de control del motor DC,
destacando especialmente el rendimiento del controlador óptimo con acción integrativa. La
incorporación del filtro de Kalman resultó fundamental para obtener estimaciones precisas de los
estados en presencia de ruido, y la fusión de sensores mediante media ponderada optimizó las
mediciones al combinar las ventajas del encoder y el tacómetro. Además, el modelado matemático
y la selección adecuada del tiempo de muestreo fueron esenciales para garantizar una respuesta
eficiente y un desempeño robusto en la simulación del sistema de lazo cerrado.
ANEXOS
Control por reubicación de polos sin AI y con AI de velocidad
clc, clear all, close all
syms s

%% Parámetros del sistema

R = 3.18; % Resistencia del devanado del motor [ohm]

L = 0.154e-3; % Inductancia devanado del motor [H]
Kt = 0.0151; % Constante de par del motor [Nm/A]
Kb = 0.0151; % Constante de fuerza contraelectromotriz [V/(rad/s)]
Bm = 0.0221e-3; % Coeficiente de fricción viscosa [Nm/(rad/s)]
J = 4.22e-7; % Inercia del motor [kg*m^2]

% Matrices de estado del sistema

a11 = -Bm/J;
a12 = Kt/J;
a21 = -Kb/L;
a22 = -R/L;

A = [a11 a12
a21 a22];

b21 = 1/L;

B = [ 0
b21];

%% Paso 1: Matriz de controlabilidad (Co)

Co = [B A*B];

if rank(Co)==2
disp('Es controlable')
else
disp('No es controlable')
end

%% Paso 2: Polinomio característico original

polos = eig(A);
poly_og = poly(polos);

a1 = poly_og(1,2);
a2 = poly_og(1,3);

%% Paso 3: Matriz W
W = [a1 1
1 0];

%% Paso 4: Matriz de Transformación (T)

T = Co*W;
%% Paso 5: Polinomio deseado
% Requerimientos
Mp = 0.05; % Factor de amortiguamiento 10%
Tst = 2; % Tiempo de estabilización (s)

% Cálculo de FACTOR DE AMORTIGUAMIENTO (e)

num = (log(Mp))^2;
den = (log(Mp))^2 + pi^2;
e = sqrt(num/den);

% Cálculo de FRECUENCIA NATURAL (wn)

wn = 4/(e*Tst);

ec = s^2 + 2*e*wn*s + wn^2;

p_des = double(solve(ec, s));

polos_des = [p_des(1,1) p_des(2,1)];

pol_des = poly(polos_des);

m1 = pol_des(1,2);
m2 = pol_des(1,3);

%% Paso 6: Matriz de ganancia (Kz)

Kz = [m2-a2 m1-a1];

%% Paso 7: Matriz de ganancia del controlador (K)

K = Kz*inv(T);

%% Paso 8: Ley de Control

% Simulación
% Tiempo de simulación
ti = 0; dt = 1e-6; tf = 10;

% Discretización
[Ak Bk] = c2d(A,B,dt);

% Condiciones iniciales
x = [0 0]';
ref = [30 0.02]';

% Inicialización de variables
k = 1;
umax = 12;
int_e = 0;

% Indicadores de desempeño
int_iae = 0;
int_ise = 0;
int_itae = 0;
int_itse = 0;

% Bucle de simulación
for tt = ti:dt:tf
x1(k) = x(1);
 x2(k) = x(2);
t(k,1) = tt;

% Error
error = ref(1)-x(1);
int_iae = int_iae + abs(error)*dt;
int_ise = int_ise + (error^2)*dt;
int_itae = int_itae + tt*abs(error)*dt;
int_itse = int_itse + tt*(error^2)*dt;

% Ley de control
u = K*(ref - x);

% Saturación de la entrada
if u > umax
u = umax;
elseif u < -umax
u = -umax;
end

V(k,1) = u;

% Actualización de estados
x = Ak*x + Bk*u;

k = k + 1;
end

IAE_RP = int_iae;
ISE_RP = int_ise;
ITAE_RP = int_itae;
ITSE_RP = int_itse;

%% Matrices de estado del sistema AMPLIADA

Ai = [ A zeros(2,1)
1 0 0 ];

Bi = [B
0];

%% Paso 1: Matriz de controlabilidad (Co)

Co = [Bi Ai*Bi Ai^2*Bi];

%% Paso 2: Polinomio característico original

polos = eig(Ai);
poly_og = poly(polos);

a1 = poly_og(1,2);
a2 = poly_og(1,3);
a3 = poly_og(1,4);

%% Paso 3: Matriz W
W = [a2 a1 1
a1 1 0
1 0 0];
%% Paso 4: Matriz de Transformación (T)
T = Co*W;

%% Paso 5: Polinomio deseado

% Cálculo de FACTOR DE AMORTIGUAMIENTO (e)

num = (log(Mp))^2;
den = (log(Mp))^2 + pi^2;
e = sqrt(num/den);

% Cálculo de FRECUENCIA NATURAL (wn)

wn = 4/(e*Tst);

p = 10*(-e*wn);

ec = (s^2 + 2*e*wn*s + wn^2)*(s-p);

p_des = double(solve(ec, s));

polos_des = [p_des(1,1) p_des(2,1) p_des(3,1)];

pol_des = poly(polos_des);

m1 = pol_des(1,2);
m2 = pol_des(1,3);
m3 = pol_des(1,4);

%% Paso 6: Matriz de ganancia (Kz)

Kz = [m3-a3 m2-a2 m1-a1];

%% Paso 7: Matriz de ganancia del controlador (K)

Ki = Kz*inv(T);

%% Paso 8: Ley de Control

% Condiciones iniciales
xi = [0 0]';
int_e = 0;

% Indicadores de desempeño
int_iae = 0;
int_ise = 0;
int_itae = 0;
int_itse = 0;

k = 1;

% Bucle de simulación
for tt = ti:dt:tf
x1i(k,1) = xi(1);
x2i(k,1) = xi(2);
t(k,1) = tt;

% Calculo del error

error = ref(1) - xi(1);
int_e = int_e + error*dt;
 int_iae = int_iae + abs(error)*dt;
int_ise = int_ise + (error^2)*dt;
int_itae = int_itae + tt*abs(error)*dt;
int_itse = int_itse + tt*(error^2)*dt;

% Paso 8: Ley de control

u = -Ki(1:2)*xi + Ki(3)*int_e;

% Saturación de la entrada
if u > umax
u = umax;
elseif u < -umax
u = -umax;
end

Vi(k,1) = u;

xi = Ak*xi + Bk*u;
k = k + 1;
end

IAE_RPi = int_iae;
ISE_RPi = int_ise;
ITAE_RPi = int_itae;
ITSE_RPi = int_itse;

figure()
subplot(311); plot(t,x1,t,x1i,'LineWidth',1.5); title('Velocidad Angular'),
xlabel('Tiempo (s)'), ylabel('w (rad/s)'), legend('Sin AI','Con AI')
subplot(312); plot(t,x2,t,x2i,'LineWidth',1.5); title('Corriente'), xlabel('Tiempo
(s)'), ylabel('I (A)'), legend('Sin AI','Con AI')
subplot(313); plot(t,V,t,Vi,'LineWidth',1.5); title('Voltaje'), xlabel('Tiempo (s)'),
ylabel('v (V)'), legend('Sin AI','Con AI')

Control por reubicación de polos sin AI y con AI de posición

clc, clear all, close all
syms s

%% Parámetros del sistema

R = 3.18; % Resistencia del devanado del motor [ohm]

L = 0.154e-3; % Inductancia devanado del motor [H]
Kt = 0.0151; % Constante de par del motor [Nm/A]
Kb = 0.0151; % Constante de fuerza contraelectromotriz [V/(rad/s)]
Bm = 0.0221e-3; % Coeficiente de fricción viscosa [Nm/(rad/s)]
J = 4.22e-7; % Inercia del motor [kg*m^2]

% Matrices de estado del sistema

a22 = -Bm/J;
a23 = Kt/J;
a32 = -Kb/L;
a33 = -R/L;

A = [0 1 0
0 a22 a23
0 a32 a33];

b31 = 1/L;

B = [ 0
0
b31];

%% Paso 1: Matriz de controlabilidad (Co)

Co = [B A*B A^2*B];

if rank(Co)==3
disp('Es controlable')
else
disp('No es controlable')
end

%% Paso 2: Polinomio característico original

polos = eig(A);
poly_og = poly(polos);

a1 = poly_og(1,2);
a2 = poly_og(1,3);
a3 = poly_og(1,4);

%% Paso 3: Matriz W
W = [a2 a1 1
a1 1 0
1 0 0];

%% Paso 4: Matriz de Transformación (T)

T = Co*W;

%% Paso 5: Polinomio deseado

% Requerimientos
Mp = 0.05; % Factor de amortiguamiento 10%
Tst = 2; % Tiempo de estabilización (s)

% Cálculo de FACTOR DE AMORTIGUAMIENTO (e)

num = (log(Mp))^2;
den = (log(Mp))^2 + pi^2;
e = sqrt(num/den);

% Cálculo de FRECUENCIA NATURAL (wn)

wn = 4/(e*Tst);

p = 10*(-e*wn);

ec = (s^2 + 2*e*wn*s + wn^2)*(s-p);

p_des = double(solve(ec, s));

polos_des = [p_des(1,1) p_des(2,1) p_des(3,1)];

pol_des = poly(polos_des);

m1 = pol_des(1,2);
m2 = pol_des(1,3);
m3 = pol_des(1,4);

%% Paso 6: Matriz de ganancia (Kz)

Kz = [m3-a3 m2-a2 m1-a1];

%% Paso 7: Matriz de ganancia del controlador (K)

K = Kz*inv(T);

%% Paso 8: Ley de Control

% Tiempo de simulación
ti = 0; dt = 5e-6; tf = 8;

% Discretización
[Ak Bk] = c2d(A,B,dt);

% Condiciones iniciales
x = [0 0 0]';
ref = [1 0 0]';

% Inicialización de variables
k = 1;
umax = 12;
int_e = 0;

% Indicadores de desempeño
int_iae = 0;
int_ise = 0;
int_itae = 0;
int_itse = 0;

% Bucle de simulación
for tt = ti:dt:tf
x1(k) = x(1);
x2(k) = x(2);
x3(k) = x(3);
t(k,1) = tt;

error = ref(1) - x(1);

int_iae = int_iae + abs(error)*dt;
int_ise = int_ise + (error^2)*dt;
int_itae = int_itae + tt*abs(error)*dt;
int_itse = int_itse + tt*(error^2)*dt;

% Ley de control
u = K*(ref - x);

% Saturación de la entrada
if u > umax
 u = umax;
elseif u < -umax
u = -umax;
end

V(k,1) = u;

% Actualización de estados
x = Ak*x + Bk*u;

k = k + 1;
end

IAE_RP = int_iae;
ISE_RP = int_ise;
ITAE_RP = int_itae;
ITSE_RP = int_itse;

%% Matrices de estado del sistema AMPLIADA

Ai = [ A zeros(3,1)
1 0 0 0 ];

Bi = [B
0];

%% Paso 1: Matriz de controlabilidad (Co)

Co = [Bi Ai*Bi Ai^2*Bi Ai^3*Bi];

%% Paso 2: Polinomio característico original

polos = eig(Ai);
poly_og = poly(polos);

a1 = poly_og(1,2);
a2 = poly_og(1,3);
a3 = poly_og(1,4);
a4 = poly_og(1,5);

%% Paso 3: Matriz W
W = [a3 a2 a1 1
a2 a1 1 0
a1 1 0 0
1 0 0 0];

%% Paso 4: Matriz de Transformación (T)

T = Co*W;

%% Paso 5: Polinomio deseado

% Cálculo de FACTOR DE AMORTIGUAMIENTO (e)

num = (log(Mp))^2;
den = (log(Mp))^2 + pi^2;
e = sqrt(num/den);
% Cálculo de FRECUENCIA NATURAL (wn)
wn = 4/(e*Tst);

p = 10*(-e*wn);

ec = (s^2 + 2*e*wn*s + wn^2)*(s-p)^2;

p_des = double(solve(ec, s));

polos_des = [p_des(1,1) p_des(2,1) p_des(3,1) p_des(4,1)];

pol_des = poly(polos_des);

m1 = pol_des(1,2);
m2 = pol_des(1,3);
m3 = pol_des(1,4);
m4 = pol_des(1,5);

%% Paso 6: Matriz de ganancia (Kz)

Kz = [m4-a4 m3-a3 m2-a2 m1-a1];

%% Paso 7: Matriz de ganancia del controlador (K)

Ki = Kz*inv(T);

%% Paso 8: Ley de Control

% Condiciones iniciales
xi = [0 0 0]';
int_e = 0;

% Indicadores de desempeño
int_iae = 0;
int_ise = 0;
int_itae = 0;
int_itse = 0;

k = 1;

% Bucle de simulación
for tt = ti:dt:tf
x1i(k,1) = xi(1);
x2i(k,1) = xi(2);
x3i(k,1) = xi(3);
t(k,1) = tt;

% Calculo del error

error = ref(1) - xi(1);
int_e = int_e + error*dt;

int_iae = int_iae + abs(error)*dt;

int_ise = int_ise + (error^2)*dt;
int_itae = int_itae + tt*abs(error)*dt;
int_itse = int_itse + tt*(error^2)*dt;

% Paso 8: Ley de control

u = -Ki(1:3)*xi + Ki(4)*int_e;

% Saturación de la entrada
 if u > umax
u = umax;
elseif u < -umax
u = -umax;
end

Vi(k,1) = u;

xi = Ak*xi + Bk*u;
k = k + 1;
end

IAE_RPi = int_iae;
ISE_RPi = int_ise;
ITAE_RPi = int_itae;
ITSE_RPi = int_itse;

figure()
subplot(321); plot(t,x1,t,x1i,'LineWidth',1.5); title('Posición Angular'),
xlabel('Tiempo (s)'), ylabel('theta (rad)'), legend('Sin AI','Con AI')
subplot(323); plot(t,x2,t,x2i,'LineWidth',1.5); title('Velocidad Angular'),
xlabel('Tiempo (s)'), ylabel('w (rad/s)'), legend('Sin AI','Con AI')
subplot(325); plot(t,x3,t,x3i,'LineWidth',1.5); title('Corriente'), xlabel('Tiempo
(s)'), ylabel('I (A)'), legend('Sin AI','Con AI')
subplot(322); plot(t,V,t,Vi,'LineWidth',1.5); title('Voltaje'), xlabel('Tiempo (s)'),
ylabel('v (V)'), legend('Sin AI','Con AI')

Control óptimo sin AI y con AI de velocidad

clc, clear all, close all
syms s

%% Parámetros del sistema

R = 3.18; % Resistencia del devanado del motor [ohm]

L = 0.154e-3; % Inductancia devanado del motor [H]
Kt = 0.0151; % Constante de par del motor [Nm/A]
Kb = 0.0151; % Constante de fuerza contraelectromotriz [V/(rad/s)]
Bm = 0.0221e-3; % Coeficiente de fricción viscosa [Nm/(rad/s)]
J = 4.22e-7; % Inercia del motor [kg*m^2]

% Matrices de estado del sistema

a11 = -Bm/J;
a12 = Kt/J;
a21 = -Kb/L;
a22 = -R/L;

A = [a11 a12
a21 a22];
b21 = 1/L;

B = [ 0
b21];

%% Paso 1: Matriz de controlabilidad (Co)

Co = [B A*B];

if rank(Co)==2
disp('Es controlable')
else
disp('No es controlable')
end

%% Paso 2: Matrices Q y R
q1 = 1e2;
q2 = 1;
Q = diag([q1 q2]);
R = [0.1];

%% Paso 3: Matriz P
P = are(A, B*inv(R)*B',Q);

%% Paso 4: Ganancia del controlador (K)

K = inv(R)*B'*P;

%% Paso 5: Ley de Control

% Tiempo de simulación
ti = 0; dt = 5e-6; tf = 5;

% Discretización del sistema

[Ak, Bk] = c2d(A, B, dt);

% Condiciones iniciales
xop = [0 0]';
ref = [1 0.005]';
umax = 12;
k = 1;

% Indicadores de desempeño
int_iae = 0;
int_ise = 0;
int_itae = 0;
int_itse = 0;

% Bucle de simulación
for tt = ti:dt:tf
x1op(k,1) = xop(1);
x2op(k,1) = xop(2);
t(k,1) = tt;

% Error
error = ref(1) - xop(1);
int_iae = int_iae + abs(error)*dt;
 int_ise = int_ise + (error^2)*dt;
int_itae = int_itae + tt*abs(error)*dt;
int_itse = int_itse + tt*(error^2)*dt;

% Ley de control
u = K*(ref - xop);

% Saturación de la entrada
if u > umax
u = umax;
elseif u < -umax
u = -umax;
end

Vop(k,1) = u;

xop = Ak*xop + Bk*u;

k = k + 1;

end

IAE_CO = int_iae;
ISE_CO = int_ise;
ITAE_CO = int_itae;
ITSE_CO = int_itse;

%% Matrices de estado del sistema AMPLIADA

Ai = [ A zeros(2,1)
1 0 0 ];

Bi = [B
0];

%% Paso 2: Matrices Q y R
q1 = 1e2;
q2 = 1;
q3 = 1e3;
Qi = diag([q1 q2 q3]);
Ri = [0.1];

%% Paso 3: Matriz P (Ec. Riccati)

P = are(Ai, Bi*pinv(Ri)*Bi',Qi);

%% Paso 4: Ganancia del controlador (Ki)

Ki = inv(Ri)*Bi'*P;

% Condiciones iniciales
xopi = [0 0]';
k = 1;
int_e = 0;

% Indicadores de desempeño
int_iae = 0;
int_ise = 0;
int_itae = 0;
int_itse = 0;

% Bucle de simulación
for tt = ti:dt:tf
x1opi(k,1) = xopi(1);
x2opi(k,1) = xopi(2);
t(k,1) = tt;

% Calculo del error

error = ref(1) - xopi(1);
int_e = int_e + error*dt;

int_iae = int_iae + abs(error)*dt;

int_ise = int_ise + (error^2)*dt;
int_itae = int_itae + tt*abs(error)*dt;
int_itse = int_itse + tt*(error^2)*dt;

% Paso 8: Ley de control

u = -Ki(1:2)*xopi + Ki(3)*int_e;

% Saturación de la entrada
if u > umax
u = umax;
elseif u < -umax
u = -umax;
end

Vopi(k,1) = u;

xopi = Ak*xopi + Bk*u;

k = k + 1;
end

IAE_COi = int_iae;
ISE_COi = int_ise;
ITAE_COi = int_itae;
ITSE_COi = int_itse;

figure()
subplot(311); plot(t,x1op,t,x1opi,'LineWidth',1.5); title('Velocidad Angular'),
xlabel('Tiempo (s)'), ylabel('w (rad/s)'),legend('CO Sin AI','CO Con AI')
subplot(312); plot(t,x2op,t,x2opi,'LineWidth',1.5); title('Corriente'),
xlabel('Tiempo (s)'), ylabel('I (A)'),legend('CO Sin AI','CO Con AI')
subplot(313); plot(t,Vop,t,Vopi,'LineWidth',1.5); title('Voltaje'), xlabel('Tiempo
(s)'), ylabel('v (V)'),legend('CO Sin AI','CO Con AI')

Control óptimo sin AI y con AI de posicion

clc, clear all, close all
syms s
%% Parámetros del sistema

R = 3.18; % Resistencia del devanado del motor [ohm]

L = 0.154e-3; % Inductancia devanado del motor [H]
Kt = 0.0151; % Constante de par del motor [Nm/A]
Kb = 0.0151; % Constante de fuerza contraelectromotriz [V/(rad/s)]
Bm = 0.0221e-3; % Coeficiente de fricción viscosa [Nm/(rad/s)]
J = 4.22e-7; % Inercia del motor [kg*m^2]

% Matrices de estado del sistema

a22 = -Bm/J;
a23 = Kt/J;
a32 = -Kb/L;
a33 = -R/L;

A = [0 1 0
0 a22 a23
0 a32 a33];

b31 = 1/L;

B = [ 0
0
b31];

%% Paso 1: Matriz de controlabilidad (Co)

Co = [B A*B A^2*B];

if rank(Co)==3
disp('Es controlable')
else
disp('No es controlable')
end

%% Paso 2: Matrices Q y R
q1 = 1e2;
q2 = 1;
q3 = 1;
Q = diag([q1 q2 q3]);
R = [1];

%% Paso 3: Matriz P
P = are(A, B*inv(R)*B',Q);

%% Paso 4: Ganancia del controlador (K)

K = inv(R)*B'*P;

%% Paso 5: Ley de Control

% Tiempo de simulación
ti = 0; dt = 5e-6; tf = 5;
% Discretización del sistema
[Ak, Bk] = c2d(A, B, dt);

% Condiciones iniciales
xop = [0 0 0]';
ref = [1 0 0]';
umax = 12;
k = 1;

% Indicadores de desempeño
int_iae = 0;
int_ise = 0;
int_itae = 0;
int_itse = 0;

% Bucle de simulación
for tt = ti:dt:tf
x1op(k,1) = xop(1);
x2op(k,1) = xop(2);
x3op(k,1) = xop(3);
t(k,1) = tt;

% Error
error = ref(1) - xop(1);
int_iae = int_iae + abs(error)*dt;
int_ise = int_ise + (error^2)*dt;
int_itae = int_itae + tt*abs(error)*dt;
int_itse = int_itse + tt*(error^2)*dt;

% Ley de control
u = K*(ref - xop);

% Saturación de la entrada
if u > umax
u = umax;
elseif u < -umax
u = -umax;
end

Vop(k,1) = u;

xop = Ak*xop + Bk*u;

k = k + 1;

end

IAE_CO = int_iae;
ISE_CO = int_ise;
ITAE_CO = int_itae;
ITSE_CO = int_itse;

%% Matrices de estado del sistema AMPLIADA

Ai = [ A zeros(3,1)
1 0 0 0 ];
Bi = [B
0];

%% Paso 2: Matrices Q y R
q1 = 1e2;
q2 = 1;
q3 = 1;
q4 = 1e3;
Qi = diag([q1 q2 q3 q4]);
R = [1];

%% Paso 3: Matriz P (Ec. Riccati)

P = are(Ai, Bi*pinv(R)*Bi',Qi);

%% Paso 4: Ganancia del controlador (Ki)

Ki = inv(R)*Bi'*P;
% Condiciones iniciales
xopi = [0 0 0]';
k = 1;
int_e = 0;

% Indicadores de desempeño
int_iae = 0;
int_ise = 0;
int_itae = 0;
int_itse = 0;

% Bucle de simulación
for tt = ti:dt:tf
x1opi(k,1) = xopi(1);
x2opi(k,1) = xopi(2);
x3opi(k,1) = xopi(3);
t(k,1) = tt;

% Calculo del error

error = ref(1) - xopi(1);
int_e = int_e + error*dt;

int_iae = int_iae + abs(error)*dt;

int_ise = int_ise + (error^2)*dt;
int_itae = int_itae + tt*abs(error)*dt;
int_itse = int_itse + tt*(error^2)*dt;

% Paso 8: Ley de control

u = -Ki(1:3)*xopi + Ki(4)*int_e;

% Saturación de la entrada
if u > umax
u = umax;
elseif u < -umax
u = -umax;
end
 Vopi(k,1) = u;

xopi = Ak*xopi + Bk*u;

k = k + 1;
end

IAE_COi = int_iae;
ISE_COi = int_ise;
ITAE_COi = int_itae;
ITSE_COi = int_itse;

figure()
subplot(321); plot(t,x1op,t,x1opi,'LineWidth',1.5); title('Posición Angular'),
xlabel('Tiempo (s)'), ylabel('theta (rad)'),legend('CO Sin AI','CO Con AI')
subplot(323); plot(t,x2op,t,x2opi,'LineWidth',1.5); title('Velocidad Angular'),
xlabel('Tiempo (s)'), ylabel('w (rad/s)'),legend('CO Sin AI','CO Con AI')
subplot(325); plot(t,x3op,t,x3opi,'LineWidth',1.5); title('Corriente'),
xlabel('Tiempo (s)'), ylabel('I (A)'),legend('CO Sin AI','CO Con AI')
subplot(322); plot(t,Vop,t,Vopi,'LineWidth',1.5); title('Voltaje'), xlabel('Tiempo
(s)'), ylabel('v (V)'),legend('CO Sin AI','CO Con AI')

Observadores por filtro de Kalman

clc, clear all, close all

% Matrices del sistema

AP = [0 1 0;
0 -0.1 0.01;
0 -0.02 -0.5];
BP = [0; 0; 0.1];
CP = [1 0 0];
DP = 0;

% Tiempo de simulación
ti = 0; tf = 10;
Tst = 1; Nr = 50;
Tau_la = -1/min(eig(AP));
T = Tau_la/8; T = Tst/Nr;

% Observador (ganancia de Kalman)

e_obs = 0.6; Tst_e = Tst / 5;
wn_obs = 4 / (e_obs * Tst_e);
polos_de_obs = [-e_obs*wn_obs+j*wn_obs*sqrt(1-e_obs^2), -e_obs*wn_obs-
j*wn_obs*sqrt(1-e_obs^2), -10*e_obs*wn_obs];
poly_de_obs = poly(polos_de_obs);

Te = [poly_de_obs(3) poly_de_obs(2) 1;
poly_de_obs(2) poly_de_obs(3) 0;
poly_de_obs(3) poly_de_obs(2) 0];
Kez = [poly_de_obs(3); poly_de_obs(2); poly_de_obs(1)] - eig(AP);
Ke = inv(Te) * Kez;

% Discretización del sistema

[Ak, Bk] = c2d(AP, BP, T);
[Aok, Bok] = c2d(AP - Ke * CP, BP, T);
[Aok, Eok] = c2d(AP - Ke * CP, Ke, T);

% Condiciones iniciales
xi = [0; 0; 0];
xoi = [0; 0; 0];
x = xi;
xo = xoi;
k = 1;

% Ruido y perturbaciones
dt = T;
t = ti:dt:tf;
nt = length(t);

rb = randn(nt, 1); % Ruido blanco - sensor

wb = randn(nt, 1); % Ruido blanco - perturbación

% Filtro pasa-bajos
wc = 1;
n = 0;
wg = 0;
for k = 1:nt
rg(k, 1) = n;
w(k, 1) = wg;
n = (1-wc*dt)*n + wc*dt*rb(k, 1);
wg = (1-wc*dt)*wg + wc*dt*wb(k, 1);
end
rg = (rg - mean(rg)) * 0.1 / std(rg);
w = (w - mean(w)) * 1 / std(w) + 1;
% Inicialización del filtro de Kalman
Pk = eye(size(AP));

% Bucle de simulación
for k = 1:nt
% Variables reales y estimadas
x1(k, 1) = x(1);
x2(k, 1) = x(2);
x3(k, 1) = x(3);
xo1(k, 1) = xo(1);
xo2(k, 1) = xo(2);
xo3(k, 1) = xo(3);
y = CP * x + rg(k, 1); % Salida medida (ruido)
% Observador con filtro de Kalman
% Predicción
xp = Ak * xo + Bk * 0 + w(k, 1);
% Corrección
Mk = Ak * Pk * Ak' + 0.1 * eye(size(AP));
Pk = inv(inv(Mk) + CP' * inv(0.01) * CP);
xo = xp + Pk * CP' * inv(0.01) * (y - CP * xp);
% Evolución del sistema
x = Ak * x + Bk * 0 + w(k, 1);
end

% Gráfica
figure
subplot(3, 1, 1);
plot(t, x1, 'b', t, xo1, 'r--');
title('Estado 1 (x1)'); grid minor;
legend('Señal original', 'Observador con KF');

subplot(3, 1, 2);
plot(t, x2, 'b', t, xo2, 'r--');
title('Estado 2 (x2)'); grid minor;
legend('Señal original', 'Observador con KF');

subplot(3, 1, 3);
plot(t, x3, 'b', t, xo3, 'r--');
title('Estado 3 (x3)'); grid minor;
legend('Señal original', 'Observador con KF');
Fusión de sensores
close all; clear; clc
LA = 0.154*10^(-3);
RA = 3.18;
KT= 15.1*10^(-3);
KB= 1/(111*(pi)/30);
J= 3.07*10^(-7);
TN=12.4*10^(-3);
WN= 5870*(pi)/30;
B=0.01*TN/WN;
% Matrices de espacio de estados
A = [0 1 0
0 -B/J KT/J
0 -KB/LA -RA/LA];
B = [0; 0; 1/LA];
C = [1 0 0];
D = 0;
% Crear el sistema de espacio de estados en MATLAB
motor_dc = ss(A, B, C, D);
% Configuración de la simulación
t = 0:0.01:5; % Tiempo de simulación de 0 a 5 segundos
V_aplicado = 12 * ones(size(t)); % Voltaje constante de 12V
% Simulación de la respuesta a un escalón de voltaje
[y, t, x] = lsim(motor_dc, V_aplicado, t);
% Gráfica de la respuesta en velocidad
figure;
plot(t, y, 'b', 'LineWidth', 1.5);
title('Respuesta en Velocidad Angular del Motor DC (Modelo 222049)');
xlabel('Tiempo (s)');
ylabel('Velocidad Angular (rad/s)');
grid on;
% Parámetros de la simulación
tiempo = t; % Tiempo de simulación de 0 a 10 segundos, con un paso de 0.01 segundos
% Velocidad real del motor (usamos una señal de ejemplo)
velocidad_real = y; %100 * sin(0.5 * tiempo); % Velocidad angular real en función del
tiempo
% Ruido de los sensores (supongamos ruido gaussiano)
ruido_sensor1 = 5 * randn(size(tiempo)); % Ruido del encoder (Sensor 1)
ruido_sensor2 = 10 * randn(size(tiempo)); % Ruido del tacómetro (Sensor 2)
% Mediciones de los sensores (velocidad real + ruido)
X1 = velocidad_real + ruido_sensor1; % Medición del encoder con ruido
X2 = velocidad_real + ruido_sensor2; % Medición del tacómetro con ruido
% Pesos para la media ponderada (basados en precisión de los sensores)
w1 = 1; % Peso asignado al Sensor 1 (encoder)
w2 = 1; % Peso asignado al Sensor 2 (tacómetro)
% Estimación combinada utilizando media ponderada
X_estimada = w1 * X1 + w2 * X2;
% Gráficas utilizando subplot para visualizar cada curva por separado
figure;
% Gráfica 1: Velocidad Real
subplot(4, 1, 1); % Dividir en una matriz 4x1 y seleccionar la primera posición
plot(tiempo, velocidad_real, 'k', 'LineWidth', 2);
title('Velocidad Real');
xlabel('Tiempo (s)');
ylabel('Velocidad (rad/s)');
grid on;
% Gráfica 2: Mediciones del Encoder (Sensor 1)
subplot(4, 1, 2); % Segunda posición
plot(tiempo, X1, 'r--');
title('Medición del Encoder (Sensor 1)');
xlabel('Tiempo (s)');
ylabel('Velocidad (rad/s)');
grid on;
% Gráfica 3: Mediciones del Tacómetro (Sensor 2)
subplot(4, 1, 3); % Tercera posición
plot(tiempo, X2, 'b--');
title('Medición del Tacómetro (Sensor 2)');
xlabel('Tiempo (s)');
ylabel('Velocidad (rad/s)');
grid on;
% Gráfica 4: Estimación Combinada mediante Media Ponderada
subplot(4, 1, 4); % Cuarta posición
plot(tiempo, X_estimada, 'g', 'LineWidth', 1.5);
title('Estimación Combinada (Media Ponderada)');
xlabel('Tiempo (s)');
ylabel('Velocidad (rad/s)');
grid on;
% Ajustar el espacio entre subplots para mejor visualización
sgtitle('Fusión de Sensores por Media Ponderada'); % Título general para todas las
subplots

Integración controlador – observador

clc, clear all, close alL
%% Datos del motor
LA = 0.154*10^(-3);
RA = 3.18;
KT = 15.1*10^(-3);
KB = 1/(111*(pi)/30);
J = 3.07*10^(-7);
TN = 12.4*10^(-3);
WN = 5870*(pi)/30;
B = 0.01 * TN / WN;
%% Modelo de estados
AP = [0 1 0
0 -B/J KT/J
0 -KB/LA -RA/LA];

BP = [0; 0; 1/LA];
CP = [1 0 0]; % Medimos únicamente la posición (x1)
DP = 0;

%% Configuración del filtro de Kalman

% Parámetros de ruido
Ruido_proceso = 1e-6 * eye(3); % Covarianza del ruido de proceso
Ruido_medida = 1e-6; % Varianza del ruido de la medida (muy pequeño para simular
ausencia de ruido)

% Matrices iniciales
Pk = eye(3); % Matriz inicial de covarianza del error
x_est = [0; 0; 0]; % Estado inicial estimado

% Discretización del sistema

Tst = 383.7e-6; Nr = 25;
T = Tst / Nr;
[Ak, Bk] = c2d(AP, BP, T);

%% Inicialización de la simulación
ti = 0; tf = 6;
dt = T; % Paso de simulación
t = ti:dt:tf;

% Condiciones iniciales
x_real = [0; 0; 0]; % Estados reales
x_integrador = 0; % Estado del integrador
u = 0; % Control inicial
r = 1; % Referencia deseada

% Variables para almacenar resultados

x1_real = []; x1_est = [];
control = []; salida = [];

%% Configuración del controlador LQR con integrador

Ai = [AP zeros(3,1); -CP 0];
Bi = [BP; 0];
Wi = [zeros(3,1); 1]; % Para referencia constante
Q = diag([1000 10 1 500]); % Penalización de estados
R = 1; % Penalización en el control
Km = lqr(Ai, Bi, Q, R);

K = Km(1, 1:3);
Ki = -Km(1, 4);

%% Simulación del sistema con observador de Kalman

for k = 1:length(t)
% Salida medida (sin ruido)
y_real = CP * x_real;
 % Predicción del filtro de Kalman
x_pred = Ak * x_est + Bk * u;
P_pred = Ak * Pk * Ak' + Ruido_proceso;

% Corrección del filtro de Kalman

Kk = P_pred * CP' / (CP * P_pred * CP' + Ruido_medida); % Ganancia de Kalman
x_est = x_pred + Kk * (y_real - CP * x_pred);
Pk = (eye(3) - Kk * CP) * P_pred;

% Controlador (basado en la posición estimada)

u = -K * x_est + Ki * x_integrador;

% Evolución del sistema real

x_real = Ak * x_real + Bk * u;

% Actualización del integrador

x_integrador = x_integrador + (r - x_est(1)) * dt; % Usando la posición estimada

% Almacenamiento de resultados
x1_real = [x1_real; x_real(1)];
x1_est = [x1_est; x_est(1)];
control = [control; u];
salida = [salida; x_est(1)];
end

%% Graficar resultados
figure;
plot(t, x1_real, 'b', t, x1_est, 'r--', 'LineWidth', 1.5);
title('Comparación: Posición Real vs Estimada'); grid on;
xlabel('Tiempo [s]'); ylabel('Posición [rad]');
legend('Posición real', 'Posición estimada');

figure;
plot(t, salida, 'b', t, r * ones(size(t)), 'r--', 'LineWidth', 1.5);
title('Posición Controlada (Estimación) vs Referencia'); grid on;
xlabel('Tiempo [s]'); ylabel('Posición [rad]');
legend('Salida estimada', 'Referencia');