Desarrollo de solvers para MPC en el lenguaje deprogramacion Julia

•

SIN SIGLA

luz angela Diaz

19/1/2024

¡Este material tiene más páginas!

Vista previa del material en texto

Equation Chapter 1 Section 1
Trabajo Fin de Grado
Ingeniería Electrónica, Robótica y Mecatrónica

Desarrollo de solvers para MPC en el lenguaje de
programación Julia
Autor: Rubén León Fuentes
Tutores: Pablo Krupa García
Ignacio Alvarado Aldea

Dpto. de Ingeniería de Sistemas y Automática
Escuela Técnica Superior de Ingeniería
Universidad de Sevilla
Sevilla, 2021

iii

Trabajo Fin de Grado
Ingeniería Electrónica, Robótica y Mecatrónica

Desarrollo de solvers para MPC en el lenguaje de
programación Julia

Autor:
Rubén León Fuentes

Tutores:
Pablo Krupa García
Ignacio Alvarado Aldea

Dpto. de Ingeniería de Sistemas y Automática
Escuela Técnica Superior de Ingeniería
Universidad de Sevilla
Sevilla, 2021

Trabajo Fin de Grado: Desarrollo de solvers para MPC en el lenguaje de programación Julia

Autor: Rubén León Fuentes
Tutores: Pablo Krupa García
Ignacio Alvarado Aldea

El tribunal nombrado para juzgar el proyecto arriba indicado, compuesto por los siguientes miembros:
Presidente:

Vocales:

Secretario:

Acuerdan otorgarle la calificación de:

El Secretario del Tribunal

Sevilla, 2021

vii

AGRADECIMIENTOS
A mis compañeros, los que estaban al principio y los que encontré por el camino. A Marina, por el apoyo y la
alegría que me ayudan a avanzar.
A mi familia, que me ha traído hasta aquí. A Ana, por ser una fuente inmensa de cariño y compañía durante
tantos y tantos días.
A Pablo, sin cuya inestimable ayuda este proyecto no sería ni la mitad de lo que es.
Todo lo que jamás se ha escrito lo ha escrito un único autor.

RESUMEN
El control predictivo por modelo (MPC) representa una de las técnicas avanzadas de control que vienen
abriéndose paso en el sector desde hace unas décadas. En este proyecto, se desarrolla un solver para la resolución
del problema matemático que surge de la implementación de MPC. El solver se desarrolla en el lenguaje
novedoso y en alza Julia y utilizada un método de optimización bien conocido en la disciplina como es ADMM.
Una vez se ha expuesto el marco teórico y se ha adecuado la formulación al problema que se desea resolver, se
desarrolla una implementación del mismo tanto en Matlab como en Julia. Posteriormente, se valida el
funcionamiento del solver a través de una serie de simulaciones frente a diferentes tipos de sistemas o parámetros
de configuración.

ABSTRACT
Since last decades, some advanced control techniques are breaking through the control systems field. Model
predictive control (MPC) represents one of them. In this project, a solver for the mathematical problem arising
from MPC implementation is developed. The solver is written in Julia, a new language currently on the rise, and
it involves the use of the optimization method ADMM.
After discussing the theoretical framework and adjusting the formulation to the problem at hand, the solver
implementation is developed both in Matlab and Julia. Afterwards, the solver performance is validated through
a series of simulations with different settings for the controlled systems and system parameters.

ÍNDICE
Agradecimientos vii
Resumen viii
Abstract ix
Índice x
Índice de Tablas xii
Índice de Ilustraciones xiii
Lista de Abreviaturas xiv
Notación xv
1 Introducción 1
1.1 Estado del arte 1
1.2 Motivación del proyecto 2
1.3 Objetivos del proyecto 2
1.4 Resumen por capítulos 3
2 Marco teórico 5
2.1 Planteamiento del problema 5
2.2 Control predictivo basado en modelo 6
2.3 Método de la alternancia de direcciones de los multiplicadores de Lagrange 8
2.4 Álgebra sparsa 9
3 Solver para MPC basado en ADMM 13
3.1 Descripción de MPC como QP 13
3.1.1 Función de coste y modelo de predicción 13
3.1.2 Desarrollo de la función de coste 15
3.1.3 Identificación de términos 17
3.2 Uso de ADMM para resolver MPC 19
3.2.1 Adaptación de la función objetivo 19
3.2.2 Desarrollo del algoritmo 20
3.3 Álgebra sparsa aplicada a ADMM 22
3.4 Implementación de un controlador MPC en bucle cerrado 23
4 Desarrollo del solver 25
4.1 Desarrollo en Matlab 25
4.1.1 Descripción del solver 25
4.1.2 Generatriz 26
4.1.3 Solver 27
4.1.4 Funciones sparsa 29
4.2 Desarrollo en Julia 32
4.2.1 Estructura de la solución propuesta en Julia 33
4.2.2 Particularidades de Julia 34
5 Resultados númericos 35
5.1 «Conic operator splitting method for convex conic problems» 35
xi

5.2 Simulaciones frente a quadprog 36
5.3 Sistema de los tanques de agua 40
5.3.1 Simulación con R = I2 41
5.3.2 Simulación con R = 10 · I2 43
5.4 Sistema de las masas conectadas por muelles 44
5.4.1 Simulación cuando no se alcanzan las restricciones y N = 10 45
5.4.2 Simulación cuando no se alcanzan las restricciones y N = 100 47
5.4.3 Simulación con restricciones activas 48
5.5 Datos estadísticos de los tiempos e iteraciones 50
6 Conclusiones 53
Anexo A: Códigos de la implementación del solver en Matlab 55
Generatriz (gen_ADMM_MPC.m) 55
Solver (solver_ADMM_MPC.m) 57
Descomposición LDL sparsa (sparsa_decomp.m) 59
Multiplicación matriz-vector sparsa (sparsa_prod.m) 59
Resolución de sistemas sparsa (sparsa_ecsys.m) 60
Coversor a CSC (full2CSC.m) 60
Coversor a CSR (full2CSR.m) 61
Simulación del sistema (sim_system.m) 62
Programa completo de la simulación (Simulaciones.m) 63
Anexo B: Códigos de la implementación del solver en Julia 65
Generación de espacio de trabajo en Julia (Generate_Julia_Workspace.m) 65
Generatriz para COSMO (gen_COSMO_MPC.m) 67
Solver ADMM (solver_ADMM.jl) 69
Funciones auxiliares para el solver ADMM (solver_ADMM.jl) 71
Simulación del sistema con el solver ADMM (solver_sim_ADMM.jl) 72
Funciones auxiliares para el solver COSMO (solver_COSMO_aux.jl) 73
Simulación del sistema con el solver COSMO (solver_sim_COSMO.jl) 74
Simulaciones para validación del solver (Tiempos.jl) 76
Referencias 79

ÍNDICE DE TABLAS

Tabla 1. Pseudocódigo de simulación 23
Tabla 2. Pseudocódigo de la función generatriz 27
Tabla 3. Pseudocódigo del solver 29
Tabla 4. Pseudocódigo de conversor a CSR 30
Tabla 5. Pseudocódigo de conversor a CSC 30
Tabla 6. Pseudocódigo de la multiplicación matriz-vector 31
Tabla 7. Pseudocódigo de la descomposición LDL 31
Tabla 8. Pseudocódigo del algoritmo QDLDL 32
Tabla 9. Datos estadísticos de los tiempos e iteraciones 51

xiii

ÍNDICE DE ILUSTRACIONES

Ilustración 1. Diagrama de sistema de control genérico basado en MPC. 7
Ilustración 2: Diagrama de control en bucle cerrado para MPC 23
Ilustración 3. Estructura de la solución propuesta en Julia 33
Ilustración 4. Estado del primer tanque de agua 37
Ilustración 5. Señal de control 1 (Grifo 1) 37
Ilustración 6. Diferencia entre la 𝑧∗calculada por quadprog y el solver ADMM propio para el sistema de los
tanques 38
Ilustración 7. Estado de la segunda masa 38
Ilustración 8. Señal de control 2 (Fuerza 2) 39
Ilustración 9. Diferencia entre la 𝑧∗calculada por quadprog y el solver ADMM propio para el sistema de los
muelles 39
Ilustración 10. Diagrama esquemático del sistema de los cuatro tanques [23] 40
Ilustración 11. Segundo estado para el sistema de los 4 tanques 41
Ilustración 12. Segunda señal de control para el sistema de los 4 tanques 42
Ilustración 13. Iteracionesdel solver para el sistema de los 4 tanques 42
Ilustración 14. Tiempos de ejecución para el sistema de los 4 tanques 43
Ilustración 15. Iteraciones del solver para el sistema de los 4 tanques (Rx10) 43
Ilustración 16. Tiempos de ejecución para el sistema de los 4 tanques (Rx10) 44
Ilustración 17. Diagrama esquemático del sistema de las masas y los muelles [24] 44
Ilustración 18. Segundo estado para el sistema de las masas y los muelles sin restricciones 46
Ilustración 19. Segunda señal de control para el sistema de las masas y los mulles sin restricciones 46
Ilustración 20. Iteraciones del solver para el sistema de las masas y los muelles sin restricciones (N=10) 47
Ilustración 21. Tiempos de ejecución para el sistema de las masas y los muelles sin restricciones (N=10) 47
Ilustración 22. Iteraciones del solver para el sistema de las masas y los muelles sin restricciones (N=100)
48
Ilustración 23. Tiempos de ejecución para el sistema de las masas y los muelles sin restricciones (N=100)
48
Ilustración 24. Segundo estado para el sistema de las masas y los muelles con restricciones 49
Ilustración 25. Segunda señal de control para el sistema de las masas y los mulles con restricciones 49
Ilustración 26. Iteraciones del solver para el sistema de las masas y los muelles con restricciones 50
Ilustración 27. Tiempos de ejecución para el sistema de las masas y los muelles con restricciones 50

file:///C:/Users/RLeon/Dropbox/Memoria/LeonFuentesRuben_def2.docx%23_Toc76401864
file:///C:/Users/RLeon/Dropbox/Memoria/LeonFuentesRuben_def2.docx%23_Toc76401871

LISTA DE ABREVIATURAS
ADMM Alternating direction method of multipliers o Método de la alternacia de
direcciones de los multiplicadores de Lagrange.

COSMO Conic operator splitting method for convex conic problemas.
Solver para optimización convexa desarrollado en Julia.

CSC Compressed sparse columna.
Formato de notación sparsa.

CSR Compressed sparse row.
Formato de notación sparsa.

LDL Nombre de un tipo de descomposición de matrices.

LQR Linear-quadratic regulator o Regulador lineal-cuadrático.

MAT Librería de Julia para la manipulación de ficheros de Matlab.

MPC Model predictive control o Control Predictivo basado en Modelo.

QDLDL Librería de Julia con herramientas para la descomposición LDL.

QP Quadratic programming o Programación cuadrática.

sparsa Dispersa.
Se dice de aquellas matrices con gran cantidad de elementos nulos.

SPCIES Toolbox de Matlab para la generación de solvers para MPC

NOTACIÓN
‖𝑢‖2 Norma euclídea del vector 𝑢

‖𝑢‖𝑄
2 Norma euclídea ponderada del vector 𝑢
‖𝑢‖𝑄
2 = 𝑢𝑇 · 𝑄 · 𝑢
𝑄 ≔ Matriz de ponderación

‖𝑢‖∞ Norma infinito del vector 𝑢
‖𝑢‖∞ = max⁡|𝑢𝑖|

𝑥𝑗|𝑘 Valor del estado 𝑥 para el instante de predicción 𝑗, en el marco de las predicciones
calculadas en el instante real 𝑘

〈𝑢, 𝑣〉 Producto escalar de los vectores 𝑢 y 𝑣

𝑥∗ Valor óptimo de la variable 𝑥
𝑓(𝑥∗) = min
𝑥
𝑓(𝑥)

1 INTRODUCCIÓN

El control predictivo basado en modelo (MPC) es un conjunto de técnicas de control basadas en el uso de un
modelo para realizar predicciones del estado futuro del sistema a controlar como parte del proceso de control
([1], [2]). Supone el marco en el que se desarrolla este proyecto y es el punto central sobre el que se mueven
todas las demás ideas planteadas. Como comienzo de este proyecto se introduce al lector en qué situación se
encuentra actualmente la tecnología de MPC en el sector de los sistemas de control. De aquí se deducen las
motivaciones y objetivos que se presuponen como guía de las tareas realizadas a lo largo del proyecto. También
se incluye un pequeño resumen por capítulos de este documento.
1.1 Estado del arte
La aplicación de sistemas de control orientados a la predicción con modelo se comienza a dar en la década de
1970, aunque algunos de los conceptos fundamentales de MPC surgen con anterioridad aplicados a otras técnicas
de control. Estos sistemas se abren paso en la industria petroquímica [3]. Utilizan modelos de respuesta a impulso
y salto escalonado y ya incluyen la implementación de restricciones en sistemas multivariables, sin embargo, no
cuentan con una teoría desarrollada que los sustente y los provea de resultados en cuanto a estabilidad o robustez
[1].
De manera independiente surgen controladores en basados modelos de función de transferencia en sistemas
monovariable que aplican ideas de control adaptativo. Son varios los diseños de controlador que aparecen bajo
esta premisa [1]. También aparecen bajo este contexto controladores basados en modelo de espacio de estados,
que permiten fácilmente ampliar a sistemas multivariable o analizar la estabilidad [4].
Aun así, la teoría no alcanzaba a proponer resultados generales de estabilidad para los controladores que
utilizaban horizonte de predicción finito. En la década de los noventa, con la aparición de nuevos controladores
centrados en espacio de estados y enfocados en este asunto, se le puso solución [5]. Actualmente, las nuevas
tecnologías que están surgiendo en el contexto del MPC se centran en obtener controladores robustos. Entre la
comunidad, hoy en día, las técnicas de control que se engloban bajo MPC son un estándar aceptado para sistemas
lineales y relativamente lentos. Los sistemas no lineales o sometidos a tiempos de muestreo muy bajos
representan el futuro del sector.
A pesar de la gran cantidad de controladores desarrollados e investigados en el mundo académico y del interés
que producen desde el punto de vista teórico, las ideas comentadas en el último párrafo han frenado
históricamente la penetración de estas tecnologías en la industria. Entre los sectores dónde más recurrente es su
uso destaca el ya mencionado como origen, sector petroquímico [6]. También se puede encontrar en la industria
aeroespacial o automovilística. Sectores donde la no linealidad de los sistemas se encuentra más presente
suponen un reto para esta técnica de control, aunque se han realizado progresos en los últimos años.
Este proyecto se enfoca en particular en un hito relativamente reciente en la historia de MPC: su aplicación en
sistemas embebidos. Esto representa un paradigma útil que permitiría llevar técnicas avanzadas de control y sus
ventajas al tipo de dispositivos controladores más presentes en la industria, como son, por ejemplo, los
Programmable Logic Controllers (PLC). En esta línea de investigación destaca la aparición de sistemas lineales
y no robustos, aunque también aparecen casos contrarios [7]. En ella se encuentran dos formas diferentes de
aplicar MPC: MPC online y MPC explícito. El primero se corresponde con las técnicas de control predictivo en
las que el problema de optimización que surge de MPC se resuelve en cada tiempo de muestreo del sistema [8].
Existen solvers diseñados específicamente con este propóstio en mente, lo que permite resultados muy eficientes.
Algunos ejemplos son μAO-MPC [9], HPMPC [10] o SPCIES [11]. La otra forma de aplicarlo, la explícita, se
aprovecha de que la acción de control óptima de un MPC lineal estándar responde a un mapa afín a trozos [12].

Introducción

2
Este mapa se puede calcular en un dispositivo ajeno al controlador, dando lugar a una look-up table. Esta tabla
permite obtener la acción de control sin resolver el problema de optimización. Esta variante destaca por su
velocidad, sin embargo, las dimensiones de la tabla crecen rápidamente con las dimensiones del problema. Esto
provoca que se requiera mucha memoria y tiempo para encontrar la solución, tanto que solo es aplicable a
sistemas pequeños con pocas restricciones.
1.2 Motivación del proyecto
El control predictivo basado en modelo es una técnica relativamente avanzada y eficiente de control que, a pesar
de estar ampliamente estudiadaen el ámbito académico, no tiene toda la penetración que cabría esperar en la
industria. Uno de los motivos a resaltar que provocan esta situación es la alta capacidad computacional que
requieren los dispositivos controladores utilizados para implementar esta técnica. Esto se debe a que MPC
(online) requiere resolver un problema de optimización matemática en cada tiempo de muestreo. Este factor es
crítico en aquellos sistemas con tiempos de muestreo muy bajos, pues los tiempos de computo hacen de cuello
de botella para cuan bajo pueden tomarse los tiempos de muestreo.
Este proyecto puede enmarcarse en el trabajo desarrollado por el profesor de la Universidad de Sevilla Pablo
Krupa, quien tiene a sus espaldas numerosas investigaciones centradas en esto. Krupa, que además es uno de lo
supervisores de este proyecto, ha escrito varios artículos académicos al respecto de la implementación de MPC
en sistemas de escasos recursos [13]. En ellos se centra en aplicar enfoques más eficientes y competitivos al
desarrollo de MPC. Además, ha desarrollado la herramienta SPCIES [11], que permite generar de forma
automatizada solvers que implementan MPC en lenguajes como C o Maltab.
Siendo el lenguaje de programación donde se implementa el solver una característica importante, se abre la
oportunidad de aplicar estos conocimientos a otros lenguajes que permitan un desarrollo aún mayor de la
eficiencia. En los últimos años se está asistiendo a una abundancia en la creación de nuevos lenguajes de carácter
más específico, con la intención de enfocarlos a tareas más concretas y con cualidades especiales que explotar
en determinados campos. Así nace Julia [14], un lenguaje caracterizado por sus logros en lo que a velocidad de
ejecución respecta. Ya existen solvers en Julia de propósito general para problemas de optimización convexa,
como es el caso de COSMO [15]. Estos solvers pueden adaptarse a MPC, pero no surgen con este en mente. Es
así que este proyecto pretende sumar un granito de arena a la implementación de solvers eficientes para MPC a
través de poner en juego estas ideas en un programa desarrollado en Julia.
1.3 Objetivos del proyecto
El objetivo central es obtener un solver funcional en el lenguaje de programación Julia. El funcionamiento del
solver se debe poder validar a través de simulación. El solver debe ser capaz de ejecutar el algoritmo ADMM
(Alternating direction method of multipliers) [16] para encontrar la solución al problema quadratic
programming (QP) al que se enfrenta, a raíz del desarrollo de la formulación elegida para MPC, en cada tiempo
de muesteo. Se pretende que el solver tenga un funcionamiento adecuado frente a diferentes tipos de sistemas.
Esto incluye sistemas con diferente número de estados o señales de control. También debe probarse el
funcionamiento con diferentes configuraciones de los parámetros que marcan las diferentes tecnologías en juego,
como son MPC o ADMM. De forma previa, ese mismo solver se validará en Matlab. Además de comprobar
que el solver puede alcanzar la referencia en una simulación, se comparará su comportamiento con el de otros
solvers disponibles, como por ejemplo el solver cuadrático nativo de Matlab o COSMO. Se espera una respuesta
similar en cuanto a estados y señales de control, pero se considera un reto particularmente difícil alcanzar la
calidad de estos solvers. Aun así, se celebrará todo lo que pueda alcanzarse en cuestión de velocidad.
Otra cuestión que tratar es la eficiencia. Se desea que el solver tenga un diseño modular que permita reducir al
mínimo el código necesario a ejecutar en cada tiempo de muestreo. Se aplica álgebra sparsa en diversas
operaciones para reducir los costes de memoria de los datos utilizados y para reducir el número de operaciones
a realizar. La implementación de estas cuestiones debe ser validada de modo que el solver eficiente que las aplica
tenga un funcionamiento correcto igual al de un solver que no las aplique. De la implementación en Julia se
espera que ofrezca unos tiempos competitivos, mejorando los resultados de Matlab y acercándose todo lo posible
a los de COSMO.

Estos representan los objetivos desde un punto de vista práctico. Aun así, no debe olvidarse que este proyecto
es un trabajo académico formativo, por lo que también se valoran los conocimientos obtenidos en las diferentes
disciplinas que lo conforman.
1.4 Resumen por capítulos
El presente documento se encuentra estructurado según se muestra a continuación. Presentados de forma
introductoria el estado del arte y la ubicación en el mismo de este proyecto, se continúa explicando las nociones
básicas de los principales conceptos que sostienen la teoría que permea este proyecto, todo ello en el Capítulo 2.
En el Capítulo 3 se aplican dichos conceptos y se desarrollan, conectándose unos a otros, de forma que se otorga
al lector el cuerpo teórico del solver desarrollado. Se incluyen todas las ecuaciones que conforman el proceso de
resolución del problema del MPC.
El Capítulo 4 consiste en la exposición de la estructura y el diseño del solver, tanto así como las decisiones que
florecieron de los retos que supuso la elaboración del mismo. Los resultados obtenidos de este proceso son
mostrados al lector a través de los pseudocódigos de cada una de las funciones que componen la herramienta
desarrollada.
Los resultados de las simulaciones con las que se ha validado el correcto funcionamiento del solver se encuentran
en el Capítulo 5. En él están expuestas las gráficas relativas al control de estados y señales de control. También
se ilustran las iteraciones y los tiempos de ejecución obtenidos. Los diferentes apartados representan a cada una
de las simulaciones, con una miríada de configuraciones diferentes, para que pueda comprobarse el
funcionamiento del solver ante situaciones interesantes de la práctica e inducir de ahí las diferentes
características de este. Finalmente se incluye una tabla con estadísticas de los tiempos de ejecución para cada
una de las simulaciones.
Las conclusiones obtenidas de dichas simulaciones son presentadas en el Capítulo 6. Ahí se reflexiona sobre lo
que ha supuesto desarrollar este proyecto, además de los beneficios obtenidos de la elaboración del mismo. Se
analizan los resultados del capítulo anterior y se desarrollan los éxitos y fracasos que suponen los resultados.

2 MARCO TEÓRICO

En este apartado se explican los conceptos teóricos sobre los que se fundamenta este proyecto, de la forma más
general posible. Es decir, en lo que concierne estas explicaciones, cada uno de los conceptos es independiente
de los demás. Por un lado, se comienza por plantear el problema, de lo que se parte y cuáles son los detalles que
se deben afrontar. Seguidamente, se habla de la técnica de control que se ha utilizado, como surge esta técnica,
por qué se caracteriza, cuál es su alcance y qué ventajas e inconvenientes aporta. La acción de control es
calculada por un solver. El algoritmo que implementa el solver es otra de las bases del proyecto. También se
responde a en qué consiste el algoritmo, cuáles son sus pasos, qué condiciones debe cumplir, en qué conceptos
matemáticos se basa. Por último, se introduce una de las herramientas utilizadas para mejorar la eficiencia del
mismo.
2.1 Planteamiento del problema
Los sistemas de control están, a grandes rasgos, compuestos de dos elementos principales: el sistema físico que
se desea controlar y el controlador encargado de realizar dicha tarea. El grueso del esfuerzo lo ocupa diseñar el
controlador, sin embargo, es necesario definir previamente que condiciones son las que describen el sistema con
el que se está trabajando. Todo sistema físico es gobernado por unas leyes que en mayor o menor grado pueden
expresarse matemáticamente a través de un modelo. Para este proyecto, a lo largo del desarrollo teórico y diseñodel sistema, no se ha trabajado con un sistema real concreto, sino sobre un modelo genérico bajo el que pueden
describirse diversos sistemas con una serie de características comunes.
El modelo al que se ha decidido recurrir es una descripción del sistema en espacio de estados dada por el
siguiente modelo lineal e invariante en el tiempo
𝑥(𝑘 + 1) = 𝐴𝑥(𝑘) + 𝐵𝑢(𝑘), 𝑥 ∈ ℝ𝑛, 𝑢 ∈ ℝ𝑚⁡ (1)
donde 𝑥(𝑘) y 𝑢(𝑘) son los estados y señales de control, respectivamente, para el instante 𝑘. Dependiendo del
sistema físico que represente pueden usarse unos u otros métodos para obtener el modelo, aunque es común
obtenerlo a través de la linealización de las ecuaciones que gobiernan el sistema aplicadas a un determinado
punto de operación. El modelo en la descripción en espacio de estados está formado por un par de ecuaciones,
sin embargo, como a lo largo de este proyecto solo se trabaja con el control de los estados y señales de control,
se obvia la ecuación correspondiente a la salida del sistema.
El modelo utilizado incluye unas restricciones, que pueden representar, por ejemplo, límites físicos que puede
alcanzar el sistema o zonas de operación del sistema que no se desea que sean sobrepasadas por alguna
determinada cuestión. Estas restricciones, tanto para los estados como para las señales de control, se modelan a
través de
𝑥 ≤ 𝑥 ≤ 𝑥, 𝑢 ≤ 𝑢 ≤ 𝑢 (2)
El objetivo de control del problema planteado es llevar el sistema a una referencia constante dada, asumiendo
que es un punto de operación del modelo, sin violar las restricciones. No violar las restricciones aplica tanto al
proceso del control (en ningún instante del control ni los estados ni las señales de control pueden salir del marco
que le imponen las restricciones) como a la referencia que se busca alcanzar (es decir, la referencia también debe
cumplir las restricciones). Esto último es lógico, pues sino la referencia sería inalcanzable.
Además de las ecuaciones mostradas, a lo largo de los desarrollos realizados en este proyecto se asumen una
serie de proposiciones de cara a definir más estrictamente las propiedades del modelo del sistema.
i. Es sistema (1) es controlable

Marco teórico

6
ii. Las retricciones (2) tienen un interior no vacío, es decir, la restricción inferior debe ser menos que la
restricción superior.
iii. La referencia (𝑥𝑟, 𝑢𝑟) es un punto de equilibrio del sistema, es decir, cumple con la ecuación
𝑥𝑟 = 𝐴𝑥𝑟 +𝐵𝑢𝑟 (3)
iv. La referencia (𝑥𝑟, 𝑢𝑟) cumple las restricciones (2), es decir, se tiene que
𝑥 ≤ 𝑥𝑟 ≤ 𝑥, 𝑢 ≤ 𝑢𝑟 ≤ 𝑢 (4)
2.2 Control predictivo basado en modelo
El control predictivo basado en modelo (MPC) es un conjunto de técnicas de control basadas en el uso de un
modelo para obtener la acción de control a través de la resolución de un problema de optimización matemática
([1], [2]). El modelo permite realizar predicciones del comportamiento del sistema en instantes futuros. El
problema de optimización consiste en la minimización de una función objetivo llamada función de coste.
Además de estas dos características, se incluye una tercera: el deslizamiento del horizonte de predicción. Esto
significa que, para cada instante en el que se realiza la predicción, se calculan predicciones para un número 𝑁
de instantes futuros, siendo 𝑁 conocido como horizonte de predicción. Para el siguiente instante, se vuelven a
culcular 𝑁 predicciones, por lo que se desplaza en una unidad el último instante calculado en la predicción con
respecto a la predicción anterior. En el contexto de este proyecto, en el que se trabaja con sistemas discretos, los
instantes de tiempo son los intervalos de muestreo. En las predicciones se incluye el cálculo de la acción de
control para dichos instantes. En cada tiempo de muestreo, la acción de control predicha para el instante actual
(la cual es solución del problema que se plantea resolver) es aplicada al sistema que se desea controlar.
Por tanto, lo más importante para las técnicas que siguen la metodología de MPC son el modelo del sistema y la
resolución del problema de optimización. La diferencia entre las distintas formulaciones y algoritmos de MPC
radica en la elección del tipo de modelo o de la función de coste, principalmente.
Este conjunto de características comunes entre todos los diferentes controladores que puede agruparse en la
familia de MPC hacen que la técnica de control seguida por estos controladores pueda describirse en los
siguientes pasos.
En primer lugar, para cada instante de tiempo, se calculan las salidas y estados para tantos instantes como indique
el horizonte de predicción. Estos cálculos dependen de dos factores, los estados pasados del sistema y las señales
de control futuras, que son aquellas que se desean calcular para posteriormente ser enviadas al sistema. En
segunda instancia, las señales de control futuras son calculadas a través de un problema de optimización, que,
por lo general, tiene como función objetivo una función cuadrática relacionada con la diferencia entre los estados
y señales de control con la referencia en el instante actual. Aunque para casos más simples puede encontrarse
una solución analítica, para los casos con restricciones, función objetivo cuadrática y modelo linear, que suponen
los más interesantes, pues la adicción de restricciones es un gran potencial de MPC, debe encontrarse la solución
a través de métodos numéricos. Finalmente, la señal de control calculada para el instante actual (el instante inicial
en las predicciones) es tomada como acción de control y es enviada al sistema. En el nuevo instante, en vez de
tomar el calculo que se hizo del estado para ese instante, se lee el estado actual del sistema, que siempre será
más fiable; y se repite el proceso desde el primer paso [1].
Para obtener un modelo de predicción pueden utilizarse diversos métodos, ya sean heurísticos o analíticos. Para
este proyecto se toma un modelo en espacio de estados. Su principal ventaja es que permite modelar sistemas
multivariable de manera sencilla. Se supone que el estado del sistema es accesible y no es necesario un
observador [1].
El concepto de utilizar un modelo para realizar predicciones del estado del sistema en instantes futuros es algo
que se hereda del linear quadratic regulator o LQR. Esta técnica recibe su nombre de utilizar sistemas lineales
con funciones de coste cuadráticas. Esta técnica propone un horizonte de eventos infinito [17]. Extender el
horizonte es útil por razones de optimalidad y factibilidad, y en ausencia de restricciones, permite encontrar una
solución analítica del problema. No obstante, las restricciones juegan un papel fundamental [2] en las

aplicaciones prácticas del control y muchos sistemas las requieren o se benefician de ellas. Es aquí donde entra
en juego MPC, aplicando un horizonte finito. Esta técnica hace resoluble el problema aún con la inclusión de
restricciones, utilizando métodos numéricos en vez de una solución analítica. Debe entonces encontrarse un
equilibrio entre tomar 𝑁 lo suficientemente altas para acercarse lo máximo posible a la solución de horizonte
infinito y lo suficientemente bajas para mantener un coste computacional adecuado.

Ilustración 1. Diagrama de sistema de control genérico basado en MPC. Figura 1.2 en [1]
Entre las ventajas que proporciona MPC frente a otras técnicas de control clásicas destacan las siguientes [1]:
• La implementación en sistemas multivariables es una extensión del caso monovariable, con lo que no
resulta especialmente compleja, a diferencia de otras técnicas de control, algunas de las cuales ni
siquiera aceptan una descripción multivariable.
• La adicción de restricciones es conceptualmente sencilla y puede incluirse fácilmente en el proceso de
diseño. Para sistemas lineales, incluir restricciones implica pasar de una solución analítica a una
numérica,generalmente de mayor costo computacional. En el caso del LQR, con horizonte infinito, la
adicción de restricciones hace irresoluble al sistema.
• Útil en sistemas en los que las referencias futuras son conocidas, como por ejemplo en el campo de la
robótica.
• Es una metodología abierta basada en unas características comunes predefinidas, lo que permite
extender su aplicación a multitud de campos.

Marco teórico

8
2.3 Método de la alternancia de direcciones de los multiplicadores de Lagrange
El método de la alternancia de direcciones de los multiplicadores de Lagrange (ADMM) es un método para la
resolución de problemas de optimización convexa. Este método surge a partir de dos métodos previos conocidos
como dual ascent y método de los multiplicadores aumentados de Lagrange [16]. Así es que combina algunas
de las características que hacen útiles a estos dos métodos. Por la parte de dual ascent, se tiene la
descomponibilidad y por la parte de los multiplicadores, las buenas propiedades para la convergencia.
Para aplicar el algoritmo se asume un problema de la forma
min
𝑧,𝑣
𝑓(𝑧) + 𝑔(𝑣) (5)
sujeto a
𝐴𝑧 + 𝐵𝑣 = 𝑑 (6)
donde 𝑓(𝑧) y 𝑔(𝑣) representan funciones convexas, 𝑧, 𝑣 ∈ ℝ𝑛 son las variables de decisión y 𝐴, 𝐵 ∈ ℝ𝑙𝑥𝑛, 𝑑 ∈
ℝ𝑙 definen una restricción lineal para las variables de decisión [16]. Sobre este problema se puede definir un
Lagrangiano aumentado
𝐿𝜌(𝑧, 𝑣, 𝜆) = 𝑓(𝑧) + 𝑔(𝑣) + 𝜆
𝑇(𝐴𝑧 + 𝐵𝑣 − 𝑑) +
𝜌
2
‖𝐴𝑧 + 𝐵𝑣 − 𝑑‖2
2 (7)
al igual que en el método de los multiplicadores. El Lagrangiano se compone de tres sumandos: la función
original, las restricciones de igualdad multiplicadas por los multiplicadores de Lagrange o variable dual 𝜆 [18]
y el término aumentado. Este último término consiste en la norma de la restricción multiplicada por un parámetro
de penalización 𝜌 [16].
A partir de aquí, se define el algoritmo como los siguientes pasos [16].
• Condición inicial: (𝑣0, 𝜆0) = (0,0) (8)
1) 𝑧𝑘+1 = argmin
𝑧
𝐿𝜌(𝑧, 𝑣
𝑘 , 𝜆𝑘) (9)
2) 𝑣𝑘+1 = argmin
𝑣
𝐿𝜌(𝑧
𝑘+1, 𝑣, 𝜆𝑘) (10)
3) 𝜆𝑘+1 = 𝜆𝑘 + 𝜌(𝐴𝑧𝑘+1 + 𝐵𝑣𝑘+1 − 𝑑) (11)

• Condición de salida:
‖𝐴𝑧𝑘+1 + 𝐵𝑣𝑘+1 − 𝑑‖
2
≤ 𝜀𝑝 (12)
‖𝑧𝑘+1 − 𝑧𝑘‖
2
≤ 𝜀𝑑 (13)
Estos pasos siguen la misma estructura que dual ascent y el método de los multiplicadores aumentados, con un
paso (o dos) de minimización de la variable (o las variables) y un paso de actualización de la variable dual. La
diferencia entre ADMM y los dos anteriores es que, en caso de aplicar estos dos a un problema de la forma (5)-
(6) la minimización de ambas variables 𝑧 y 𝑣 se realiza de manera simultánea, mientras que en ADMM ese paso
se desdobla en dos pasos diferentes.
Se puede observar en la práctica que las propiedades de convergencia de ADMM no lo hacen un algoritmo
especialmente rápido, pero sí lo suficiente como para alcanzar resultados de una precisión considerable en unas
decenas de iteraciones para la mayoría de aplicaciones prácticas [16].

2.4 Álgebra sparsa
El álgebra sparsa es un tipo de notación que permite reducir el espacio de almacenamiento y los tiempos de
cómputo en cálculos relacionados con matrices de gran tamaño y gran cantidad de elementos nulos. Se trata de
una herramienta matemática que está estrechamente relacionada y aporta utilidad en el cálculo numérico y en
las ciencias de la computación. Una de las premisas clave para este proyecto es la eficiencia y en ese sentido es
donde encuentra un lugar al uso de este tipo de notación, diseñada específicamente con la eficiencia como
objetivo.
Existen diversos métodos para describir matrices de forma sparsa. En general, todos consisten en guardar las
componentes no nulas y los índices que las sitúan dentro de la matriz en diversos vectores. A lo largo de este
proyecto se utilizan dos formatos: compressed sparse row (CSR) y compressed sparse column (CSC) [19]. Las
sucesivas explicaciones se realizan sobre el primer formato, dado que el segundo es idéntico, cambiando la
función realizada por las filas a las columnas. En primer lugar, debe definirse como funciona el formato y,
seguidamente, se explican los dos tipos de operaciones que se han realizado a lo largo del proyecto con este tipo
de matrices: el producto matriz-vector y la resolución de sistemas de ecuaciones de la forma 𝐴𝑥 = 𝑏. Nótese
que es un tema estrechamente relacionado con la computación, así que se debe mencionar que las explicaciones
dadas se basan en indexación con comienzo en el uno.
Para el formato CSR, toda matriz queda descrita con un conjunto de tres vectores. El primero de ellos, vector de
valores, contiene todos los valores distintos de cero, de izquierda a derecha y de arriba abajo. El segundo de
ellos, vector de columnas, contiene los índices de la columna de cada uno de esos valores. Es decir, el primer
elemento de este vector tiene almacenado la columna en la que se encuentra el primer elemento del vector de
valores, y así sucesivamente. El último vector es el que rompe esta dinámica, el vector de filas. En este vector,
la n-ésima componente del vector guarda información de la n-ésima fila. La información que guarda es en qué
componente del vector de valores comienza dicha fila. El vector de filas tiene un tamaño igual al número de filas
más uno. El último valor que queda almacenado en él es el número total de elementos no nulos más uno (o visto
de otra forma, la longitud del vector de valores más uno). Esto cobra sentido cuando se reconstruye la matriz. El
ejemplo (14) ilustra cómo funciona el formato. Los subíndices marcan la posición que ocupa el elemento en el
vector o matriz correspondiente.

𝑀 = (
21,𝟏 31,𝟐 0
0 42,𝟐 0
0 0 63,𝟑
0 0 0
0 0 0
0 72,𝟓 0
73,𝟒 83,𝟓 0
0 0 104,𝟔
)
(14)

𝑣𝑒𝑐𝑡𝑜𝑟𝑣𝑎𝑙𝑜𝑟𝑒𝑠 = (2𝟏 32 4𝟑 74 6𝟓 76 87 10𝟖)
𝑣𝑒𝑐𝑡𝑜𝑟𝑐𝑜𝑙𝑢𝑚𝑛𝑎𝑠 = (1 2 2 5 3 4 5 6)
𝑣𝑒𝑐𝑡𝑜𝑟𝑓𝑖𝑙𝑎𝑠 = (1 3 5 8 9)
El formato CSR permite leer fácilmente los elementos de la matriz por filas. Para reconstruir una matriz a partir
de su descripción sparsa, es necesario extraer los datos fila a fila y luego concatenar esas filas en una única
matriz. Por eso, este formato resulta especialmente útil para la multiplicación matriz-vector. En ese tipo de
operación, la componente n-ésima del vector resultado viene dada por el producto escalar de la n-ésima fila de
la matriz por el vector multiplicado. Es así que conviene poder acceder a la matriz por filas. Para la multiplicación
vector-matriz, resulta más favorable el formato CSC. Para reconstruir las filas se sigue el siguiente proceso. Se
selecciona la componente del vector de filas correspondiente a la fila que se desea reconstruir como comienzo y
la siguiente componente como fin. Del vector de valores, se toman las componentes de índice comprendido entre
el comienzo y el fin. Así se obtienen todas las componentes no nulas de esa fila. Se repite el proceso en el vector
de columnas y así se obtiene la posición exacta que ocupa cada valor dentro de la fila. En el ejemplo (15) se
reconstruye la fila 2 (𝐹𝑀2) de (14). Si se quiere construir la matriz, el resto de elementos se rellenan con ceros.
Si se desea hacer el producto matriz-vector, cada elemento no nulo se multiplica por el elemento correspondiente
del vector, según indique el índice de columnas.

Marco teórico

10
𝑣𝑒𝑐𝑡𝑜𝑟𝑓𝑖𝑙𝑎𝑠 = (1 𝟑2 𝟓2+1 8 9)
(15)
𝑣𝑒𝑐𝑡𝑜𝑟𝑣𝑎𝑙𝑜𝑟𝑒𝑠 = (21 32 {4𝟑 74 }6𝟓 76 87 108)
𝑣𝑒𝑐𝑡𝑜𝑟𝑐𝑜𝑙𝑢𝑚𝑛𝑎𝑠 = (11 22 {2𝟑 54 }3𝟓 46 57 68)

𝐹2𝑀 = (0 42,𝟐 0 0 72,𝟓 0)
En el ejemplo (16) puede observarse como se realiza dicha multiplicación. La matriz original de (14) se
multiplica por un vector 𝑣, dando como resultado un vector 𝑠. El segundo elemento del vector resultado es igual
a la multiplicación escalar de la segunda fila de la matriz por el vector𝑣. Los valores de la fila 2 ya se extrajeron
en (15). Cada uno de esos valores se multiplica por el elemento del vector que ocupa la posición que indique el
vector de columnas según (15). Importante destacar que las multiplicaciones por los elementos nulos no se
realizan, pues ya se sabe de antemano que no aportan nada al resultado final Concluyendo, se suman los
resultados de los productos y ese es el resultado final para el segundo elemento de la solución.
𝑀4𝑥6 · 𝑣6𝑥1 = 𝑠4𝑥1
con⁡𝑣 = (1 2 3 4 5 6)𝑇
(16)

𝑠2 = 〈𝐹2𝑀 , 𝑣〉 = (· 4 · · 7 ·) ·
(

·
𝑣𝟐
·
·
𝑣𝟓
· )

4 · 𝑣𝟐 + 7 · 𝑣𝟓 = 4 · 2 + 7 · 5 = 43

𝑀 · 𝑣 = (
2 3 0
0 4 0
0 0 6
0 0 0
0 0 0
0 7 0
7 8 0
0 0 10
)
4𝑥6
·
(

1
2
3
4
5
6)

6𝑥1
= (
·
432
·
·
)
4𝑥1

El último aporte de la notación sparsa es la resolución de sistemas de ecuaciones. Las operaciones requeridas en
este caso son más complejas y están basadas en QDLDL. Esto es una herramienta que permite la resolución de
sistemas del tipo 𝐴𝑥 = 𝑏 a través de descomposición LDL, desarrollado por el University of Oxford Control
Group como parte del proyecto del solver OSQP [20]. Dentro de la librería existen hasta cinco funciones para
realizar manipulación de sistemas LDL. Una de ellas permite resolver sistemas del tipo 𝐿𝐷𝐿𝑇𝑥 = 𝑏, y el
algoritmo que implementa es la principal influencia para la metodología aplicada en este proyecto.
Por tanto, el algoritmo utilizado para resolver los sistemas requiere de un paso intermedio. El algoritmo se aplica
sobre la descomposición LDL de la matriz de coeficientes. Esta descompisicón está relacionada con la
descomposición de Cholesky, que permite reescribir una matriz como el producto de una matriz triangular y su
traspuesta conjugada. En el caso de que los elementos de la matriz sean reales, como es este caso, la traspuesta
conjugada es igual a la traspuesta. Las matrices definidas positivas tienen por propiedad que la descomposición
de Cholesky es única. En el caso de la decomposición LDL, la matriz triangular 𝐿 tiene una diagonal unitaria,

es decir, todos sus elementos valen uno. Pero esta no es la única matriz que aparece en el producto. También
aparece una nueva matriz diagonal 𝐷 entre medias del producto de las matrices triangulares. Es decir, la matriz
𝐷 queda multiplicada por 𝐿 a su izquierda y por 𝐿 traspuesta a su derecha. De ahí el nombre de este tipo de
descomposición.
Resulta útil para cuestiones de implementación relacionar ambas descomposiciones, así puede realizarse la LDL
a partir de una Cholesky dada. Sea
𝑀 = 𝑅 · 𝑅𝑇 (17)
la descomposición de Cholesky de una matriz genérica definida positiva y 𝑆 una matriz diagonal con los
elementos de la diagonal principal de 𝑅. Sea
𝑀 = 𝐿 · 𝐷 · 𝐿𝑇 (18)
la descomposición LDL de la misma matriz. Si se desarrolla (18) de la forma indicada

𝑀 = 𝐿𝐷𝐿𝑇 = 𝐿 · 𝐷
1
2⁄ · (𝐷
1
2⁄ )
𝑇
· 𝐿𝑇 = 𝐿𝐷
1
2⁄ · (𝐷
1
2⁄ 𝐿)
𝑇
(19)
y teniendo en cuenta que al ser 𝑀 definida positiva ambas descomposiciones son únicas, se obtiene

𝑀 = 𝐿𝐷
1
2⁄ · (𝐷
1
2⁄ 𝐿)
𝑇
= 𝑅 · 𝑅𝑇 ⁡⇒ 𝑅 = 𝐿𝐷
1
2⁄ (20)
Esta conclusión permite obtener la descomposición Cholesky a partir de la LDL. Es más común que se requiera
el fenómeno contrario y, por tanto, el resultado queda mejor expresado como
𝑅 = 𝐿𝐷
1
2⁄ ⁡⇒ ⁡ { ⁡𝐷 = 𝑆
2
⁡𝐿 = 𝑅𝑆−1
(21)

3 SOLVER PARA MPC BASADO EN ADMM

Los dos pilares sobre los que se fundamenta este proyecto son los sistemas de control y la optimización
matemática, a través de las técnicas de MPC y ADMM, respectivamente. Sin embargo, podría parecer que se
trata de dos áreas de conocimiento ajenas entre sí. A lo largo de este apartado, se pone de manifiesto la forma en
la que ambas quedan relacionadas. En concreto, como puede describirse el problema surgido de MPC como un
problema de optimización, que puede ser resuelto mediante el algoritmo ADMM. A través de la descripción de
este proceso, también se aclaran el planteamiento, las particularidades y las asunciones del problema exacto que
se está resolviendo, tanto así como la formulación y desarrollo matemáticos utilizados. También se describe
como el resto de herramientas que se han utilizado facilitan este desarrollo y lo hacen más eficiente.
3.1 Descripción de MPC como QP
En el núcleo del MPC se encuentra la función de coste. Esta expresión describe las diferencias entre las variables
en el estado actual del sistema y los valores de referencia y es sobre ella que interesa actuar. El objetivo del
control no deja de ser que las acciones de control y los estados se ajusten a la referencia de interés, y siendo la
función de coste la descripción de la diferencia de ambas, el ideal sería que el resultado de esta función, es decir,
el coste, fuese nulo. No siempre esto es posible, por tanto, es deseable encontrar la actuación sobre el sistema
que minimice el coste. Es decir, el objetivo del control puede reducirse a encontrar unos valores que minimicen
una determinada expresión matemática, lo que no deja de ser un problema de optimización. Particularizando, la
función de coste resulta en una función cuadrática, por lo que el problema resultante puede ser descrito con la
formulación de quadratic programming o QP.
3.1.1 Función de coste y modelo de predicción
El conjunto de ecuaciones que representa el planteamiento del problema es (1)-(4). A partir de él se ha
seleccionado la siguiente descripción de la función de coste [17].

𝐽 = min
𝑥,𝑢
[‖𝑥𝑁 − 𝑥𝑟‖𝑃
2 +∑‖𝑥𝑗 − 𝑥𝑟‖𝑄
2
+ ‖𝑢𝑗 − 𝑢𝑟‖𝑅
2
𝑁−1
𝑗=0
] (22)
sujeto a
• Modelo de
predicción:
𝑥𝑗+1 = 𝐴𝑥𝑗 + 𝐵𝑢𝑗 (23)

• Restricciones:
𝑥 ≤ 𝑥𝑗 ≤ 𝑥,⁡⁡⁡⁡⁡⁡⁡⁡⁡𝑗 = 1…𝑁
𝑢 ≤ 𝑢𝑗 ≤ 𝑢, 𝑗 = 1…𝑁
(24)
• Condición
inicial:
𝑥0 = 𝑥(𝑘) (25)
• Restricción
final:
𝑥𝑁 ∈ Ω → 𝑥𝑁 = 𝑥𝑟 (26)
A continuación, se explican en detalle los diferentes conceptos, elementos y símbolos que componen estas
expresiones. El primer bloque sobre el que se construye esta función es la diferencia entre valores actuales y la

Solver para MPC basado en ADMM

14
referencia. La idea es que un único valor de coste represente las diferencias de todas las componentes de los
estados y señales de control con sus respectivas referencias, por tanto, la diferencia se realiza a través de la norma
euclídea ponderada
‖𝑥𝑗 − 𝑥𝑟‖𝑄
2
⁡; ‖𝑢𝑗 − 𝑢𝑟‖𝑅
2
(27)
El hecho de ponderarla permite incluir parámetros de diseño, 𝑄 y 𝑅, que regulen cuanto peso tiene cada una de
esas componentes sobre el valor del coste, en caso de que fuese necesario. Estos parámetros se llaman matriz de
ponderación y se asume que deben ser semidefinidas positivas, aunque en la mayoría de los casos se tomen
como estrictamente definidas positivas. Una condición para ello es que sean simétricas, aunque en las
simulaciones realizadas para validación se toman directamente diagonales. Esta condición es usada más adelante
para simplificar algunas expresiones. La norma se aplica de manera independiente a los vectores de estados y
señales de control y luego se suman ambas,
‖𝑥𝑗 − 𝑥𝑟‖𝑄
2
+ ‖𝑢𝑗 − 𝑢𝑟‖𝑅
2
(28)
pues el coste tiene en cuenta a ambas de manera conjunta. Sin embargo, se reserva el estado correspondiente al
horizonte de predicción para calcularlo aparte, con una matriz de ponderación diferente (𝑃), según
‖𝑥𝑁 − 𝑥𝑟‖𝑃
2 (29)
Esta decisión se toma por cuestiones de estabilidad y optimalidad. Tratándose de control predictivo, el coste no
debe tener en cuenta solo los valores en el estado actual, sino todos los valores en los instantes dentro del
horizonte de predicción 𝑁. Por ello, el desarrollo se aplica desde el instante actual hasta el horizonte de
predicción, a través del sumatorio

‖𝑥𝑁 − 𝑥𝑟‖𝑃
2 +∑‖𝑥𝑗 − 𝑥𝑟‖𝑄
2
+ ‖𝑢𝑗 − 𝑢𝑟‖𝑅
2
𝑁−1
𝑗=0
(30)
Se asumeque el horizonte de predicción tiene que ser mayor que la dimensión del vector de estados. En general,
estos son los componentes que construyen la función. En esta función, las variables de decisión, de lo que
depende el coste, son los valores de estados y señales de control en los diferentes instantes. Por tanto, esto es lo
que interesa averiguar, determinando los valores que hacen mínimo el coste. El valor del coste mínimo en sí no
es particularmente interesante. El objetivo último es obtener los valores de las señales de control para el instante
actual que hacen el coste mínimo. Aun así, primero se debe resolver el problema de optimización, con todas las
variables, y para ello se parte de la función de coste mencionada en (22)-(26).
Cabe resaltar un detalle fundamental. Puede observarse que la notación utilizada para representar el valor de
estados y señales de control en cada instante es diferente en el modelo de predicción y en el modelo del sistema.
Esto se debe a que el modelo sobre el que opera la función de coste no es el modelo del sistema, sino el modelo
de predicción. Aunque estos dos modelos tengan la misma apariencia matemática, no son iguales
conceptualmente hablando. Para un determinado instante de funcionamiento 𝑘, el sistema físico se encuentra en
un estado representado en el modelo del sistema por 𝑥(𝑘). Para ese mismo instante, durante el cálculo de control,
el controlador debe calcular los valores de estado para los siguientes 𝑁 instantes. Esos valores de los estados son
llamados 𝑥𝑗|𝑘 , 𝑥𝑗+1|𝑘…𝑥𝑗+𝑁|𝑘. Debe observarse que en la notación 𝑥𝑗|𝑘, el subíndice 𝑗 hace referencia al valor
que ocupa 𝑥 en las predicciones, mientras que el subíndice 𝑘 hace referencia al instante de funcionamiento en
el que se están realizando las predicciones. Una vez aplicada la señal de control en ese instante, el sistema
conjunto avanza hacia el instante 𝑘 + 1, con los estados 𝑥(𝑘 + 1). Nótese que 𝑥(𝑘 + 1) y 𝑥𝑗+1|𝑘 no son iguales
a pesar de hacer referencia al mismo instante de funcionamiento, pues no existe un modelo perfecto y no se
puede predecir el estado del sistema con total exactitud. En este nuevo instante se vuelve a realizar el control y
de nuevo se calculan los estados 𝑥𝑗|𝑘+1, 𝑥𝑗+1|𝑘+1…𝑥𝑗+𝑁|𝑘+1. Con la notación utilizada, el nuevo 𝑥𝑗|𝑘+1 hace
referencia al instante real 𝑘 + 1, pues el modelo de predicción se aplica de nuevo independientemente de lo que
sucediese en el instante anterior. Por describirlo mediante una analogía, es como un par de bucles anidados, en
el que el índice del segundo bucle se reinicia cada vez que el primero comienza una iteración nueva. Es decir, el

nuevo 𝑥𝑗+1|𝑘+1 hace referencia al mismo instante real (𝑘 + 2) del sistema físico que el 𝑥𝑗+2|𝑘 calculado en la
primera predicción. Sin embargo, el cálculo del nuevo 𝑥𝑗+1|𝑘+1 está basado en el valor de 𝑥(𝑘 + 1), que se
ajusta mucho más a la realidad que el antiguo valor de 𝑥𝑗+1|𝑘, en el cual se basa el cálculo del antiguo 𝑥𝑗+2|𝑘.
Teniendo en cuenta que los valores útiles para el sistema físico son los calculados para la predicción
correspondiente al instante en el que se realizan las predicciones, el cálculo de los demás puede parecer
innecesario. Sin embargo, debe recordarse del apartado 2.2 que este planteamiento viene heredado del
planteamiento del LQR y, por tanto, a mayor longitud del horizonte de predicción, mas cercano es el resultado
al de horizonte infinito. Esto implica mejor optimalidad y factibilidad.
Realizada la apreciación, debe discutirse también la relación entre el modelo de predicción y el modelo del
sistema (el modelo de predicción no deja de ser también un modelo del sistema, pero con esta nomenclatura se
hace referencia al modelo utilizado como planteamiento del problema). Ambos modelos siguen notaciones
diferentes porque conceptualmente representan cosas diferentes, sin embargo, tienen una estructura similar. Al
que se hace referencia como modelo del sistema solo es una abstración matemática utilizada para los desarrollos
teóricos que se están realizando. El modelo de predicción es una herramienta utilizada por MPC para poder
controlar el sistema. Para una aplicación concreta con un sistema físico real, ambos modelos representarían el
mismo sistema. Es más, si el sistema físico no es real ni simulado, ambos modelos son directamente iguales. El
modelo de predicción es único, pero se repite de manera independiente (de un instante a otro) para cada instante
que transcurre en el modelo del sistema. Tal y como se puede observar en el ejemplo del párrafo anterior, para
cada instante real, el valor inicial del modelo de predicción 𝑥0 se corresponde con el valor actual del modelo del
sistema 𝑥(𝑘) (25). También debe cumplirse una restricción final, que para este proyecto se toma lo más sencilla
posible. Esto quiere decir que el estado al final del horizonte de predicción 𝑥𝑁 debe pertenecer a un conjunto
convexo determinado Ω. El conjunto más sencillo se traduce a que dicho estado debe ser igual a la referencia,
conclusión lógica (26).
3.1.2 Desarrollo de la función de coste
La función de coste ha sido presentada de forma que los términos que la componen reflejen de manera clara los
conceptos relacionados con sistemas de control que la constituyen. En la expresión (22) puede entenderse
claramente que ideas hay tras las normas o el porqué se está realizando un sumatorio. Si se desarrollan las
normas, pues observarse sin complicación que se trata de una función cuadrática. No obstante, en la literatura
matemática, estas funciones suelen escribirse en otro tipo de expresión, para facilitar la resolución de las mismas.
La manera más común de representar las funciones cuadráticas es a través de expresiones matriciales en un
único vector de variables. Tal y como se encuentra la formulación, existen diferentes vectores que representan
las variables de decisión del problema. Es por tanto necesario desarrollar un poco más la expresión (22) y agrupar
todas las variables para llegar a una formulación matricial diferente del problema.
En general, en optimización, los problemas de minimización de funciones cuadráticas con restricciones lineales
responden al nombre QP. Este tipo de problemas cuenta con diversas formas de expresarlos. Para el desarrollo
de este proyecto se ha elegido la formulación

min
𝑧
1
2
𝑧𝑇𝐻𝑧 + 𝑞𝑇𝑧 (31)
sujeto a
𝐴𝑒𝑞𝑧 = 𝑏 (32)
𝐺𝑧 ≤ 𝑑 (33)
donde 𝑧 ∈ ℝZ es la variablesde decisión, 𝐻 ∈ ℝZxZ es la matriz Hessiana, 𝑞 ∈ ℝZ son los coeficientes del
término lineal, 𝐴𝑒𝑞 ∈ ℝ
𝑙𝑥𝑍, 𝑏 ∈ ℝ𝑙 definen una restricción lineal para las variables de decisión y 𝐺 ∈ ℝ𝑝𝑥𝑍, 𝑑 ∈
ℝ𝑝 definen una restricción de desigualdad para las variables de decisión. Se trata de una expresión bastante
genérica del mismo. Aun así, puede realizarse un ajuste que facilite más la caracterización del MPC. La
restricción de desigualdad (33) representa una delimitación del espacio que puede tomar la variable del
problema. Esa inecuación concreta es útil a la hora de representar un lugar geométrico genérico, incluidos
volúmenes más complejos. Sin embargo, las restricciones que se han planteado como enunciado del MPC, que

Solver para MPC basado en ADMM

16
son las que definen esta condición de desigualdad, representan un lugar geométrico más sencillo: un cuadrado,
cubo o equivalente en dimensiones mayores. Por tanto, las condiciones de desigualdad pueden mantenerse tal y
como quedan planteadas en el MPC, resultado así el problema QP en

min
𝑧
1
2
𝑧𝑇𝐻𝑧 + 𝑞𝑇𝑧 (34)
sujeto a
𝐴𝑒𝑞𝑧 = 𝑏 (35)
𝑧 ≤ 𝑧 ≤ 𝑧 (36)
Quedan así tres elementos entre los que establecer la conexión: las variables del problema, la función objetivo y
la condición de igualdad. La manera más cómoda de hacerlo será identificando términos. Para ello, debe
realizarse un paso previo, desarrollarla función (22).
En primer lugar, la norma ponderada puede ser desarrollada como
‖𝑥𝑗 − 𝑥𝑟‖𝑄
2
= (𝑥𝑗 − 𝑥𝑟)
𝑇
𝑄(𝑥𝑗 − 𝑥𝑟) = (𝑥𝑗
𝑇𝑄 − 𝑥𝑟
𝑇𝑄)(𝑥𝑗 − 𝑥𝑟)
= 𝑥𝑗
𝑇𝑄𝑥𝑗 − 𝑥𝑗
𝑇𝑄𝑥𝑟 − 𝑥𝑟
𝑇𝑄𝑥𝑗 + 𝑥𝑟
𝑇𝑄𝑥𝑟
= 𝑥𝑗
𝑇𝑄𝑥𝑗 − 2 · 𝑥𝑟
𝑇𝑄𝑥𝑗 + 𝑥𝑟
𝑇𝑄𝑥𝑟
(37)
Debe notarse que, de los tres términos resultantes, uno de ellos solo depende de la referencia y de la matriz de
ponderación. En consecuencia, este término afecta al valor final del coste, pero no a que valor de las variables
de decisión lo hace mínimo y, por tanto, podemos despreciarlo del desarrollo, dejando el resultado en dos
sumandos. Este mismo procedimiento puede aplicarse a las otras dos normas ponderadas que aparecen en (22).
La función de coste queda entonces como

𝐽 = min
𝑥,𝑢
[(𝑥𝑁
𝑇𝑃𝑥𝑁 − 2𝑥𝑟
𝑇𝑃𝑥𝑁) + ∑(𝑥𝑗
𝑇𝑄𝑥𝑗 − 2𝑥𝑟
𝑇𝑄𝑥𝑗) + (𝑢𝑗
𝑇𝑅𝑢𝑗 − 2𝑢𝑟
𝑇𝑅𝑢𝑗)
𝑁−1
𝑗=0
] (38)
Si se deshacen los paréntesis, pueden reagruparse los términos de la siguiente manera. Por un lado, se agrupan
los términos en los que aparece una matriz ponderada multiplicada a su derecha por un vector de variables y a
su izquierda por la traspuesta de ese mismo vector. Por otro lado, se agrupan los términos que tienen una variable
multiplicando a la derecha y a un vector de referencias a la izquierda. Es más, el dos que multiplica a todos los
elementos de este segundo grupo puede dejarse fuera del paréntesis como factor común. El resultado es

𝐽 = min
𝑥,𝑢
[(𝑥𝑁
𝑇𝑃𝑥𝑁 +∑ 𝑥𝑗
𝑇𝑄𝑥𝑗 + 𝑢𝑗
𝑇𝑅𝑢𝑗
𝑁−1
𝑗=0
)
− 2(𝑥𝑟
𝑇𝑃𝑥𝑁 +∑ 𝑥𝑟
𝑇𝑄𝑥𝑗 + 2𝑢𝑟
𝑇𝑅𝑢𝑗
𝑁−1
𝑗=0
)]
(39)
La forma de esta expresión es parecida a la forma de la función objetivo del problema QP, sin embargo, es
necesario reagrupar más los términos, de modo que se encuentren unidas todas las variables cuadradas por un
lado y las lineales por otro. Para ello, debe desarrollarse el sumatorio de la siguiente manera
𝐽 = min
𝑥,𝑢
[𝑥0
𝑇𝑄𝑥0 +⋯+ 𝑥𝑁−1
𝑇 𝑄𝑥𝑁−1 + 𝑥𝑁
𝑇𝑃𝑥𝑁 + 𝑢0
𝑇𝑅𝑢0 +⋯+ 𝑢𝑁−1
𝑇 𝑅𝑢𝑁−1
− 2(𝑥𝑟
𝑇𝑄𝑥0 +⋯+ 𝑥𝑟
𝑇𝑄𝑥𝑁−1 + 𝑥𝑟
𝑇𝑃𝑥𝑁 + 𝑢𝑟
𝑇𝑅𝑢0 +⋯
+ 𝑢𝑟
𝑇𝑅𝑢𝑁−1)]
(40)

3.1.3 Identificación de términos
En este punto, ya puede redefinirse la variable del problema, de modo que se ajuste a la del problema QP.
Simplemente se construye un nuevo vector que contenga a todas las variables, los estados y señales de control
para cada instante. No existe una forma única de distribuir las variables dentro del vector. Para este proyecto, se
ha seleccionado la siguiente

𝑧 =
(

𝑢0
𝑥1
𝑢1
𝑥2
𝑢2
⋮
𝑥𝑁−1
𝑢𝑁−1
𝑥𝑁 )

(41)
Esta forma hace más fácil la construcción del resto de elementos. Además, las variables que interesan desde el
punto de vista del control (las señales de control en el instante actual) ocupan las primeras posiciones del vector.
Si se reescriben todos los sumandos en (40) como una operación matricial, de modo que las variables de la
ecuación resulten en (41), las matrices de ponderación formarán una nueva matriz. Esta matriz es el equivalente
de la matriz 𝐻 del problema QP. Por tanto, la formación de esta matriz queda como

𝐻 =
(

𝑅
𝑄
𝑅
𝑄
𝑅
⋱
𝑄
𝑅
𝑃)

(42)
Para los términos lineales, puede aplicarse el mismo razonamiento, siendo el resultado

𝑞 = −
(

𝑅𝑢𝑟
𝑄𝑥𝑟
𝑅𝑢𝑟
⋮
𝑄𝑥𝑟
𝑅𝑢𝑟
𝑃𝑥𝑟)

(43)
Cabe comentar un par de detalles. En primer lugar, mientras el término lineal del MPC aparecía restando, en la
formulación QP aparece sumando, por lo que todas las componentes del vector aparecen con el signo cambiado.
En segundo lugar, al ser las matrices ponderadas simétricas por definición, es equivalente multiplicarles el vector
traspuesto de referencia por la izquierda, que multiplicarles el mismo por la derecha sin trasponer. Además, toda
la expresión puede dividirse entre dos sin alterar la solución óptima. Esto se debe a que la solución óptima que
se busca es la variable de decisión óptima, y no el valor óptimo de la función (que sí quedaría alterado).
Una vez se ha formado el vector de variables, atendiendo a (36), los vectores de cota superior e inferior resultan
en

Solver para MPC basado en ADMM

𝑧 =
(

𝑢
𝑥
⋮
𝑢
𝑥)

𝑧 =
(

𝑢
𝑥
⋮
𝑢
𝑥)

(44)
El único elemento restante para terminar de adaptar MPC a QP son las restricciones de igualdad.
En el problema planteado, las restricciones de igualdad vienen dadas por la condición inicial del problema (25).
Si se sustituye esta expresión en el modelo de predicción se obtiene
𝑥1 = 𝐴𝑥0 + 𝐵𝑢0 ⁡→ ⁡ 𝑥1 = 𝐴𝑥(𝑘) + 𝐵𝑢0 (45)
Redistribuyendo los miembros de la ecuación, con los datos numéricos en el primer miembro y las variables en
el segundo
−𝐴𝑥(𝑘) = 𝐵𝑢0 − 𝑥1 (46)
Esta ecuación puede construirse para todos los instantes del modelo de predicción hasta el horizonte de eventos.
Para el segundo instante, la ecuación es
𝑥2 = 𝐴𝑥1 + 𝐵𝑢1 ⁡→ ⁡𝐴𝑥1 + 𝐵𝑢1 − 𝑥2 = 0 (47)
Este proceso puede aplicarse a todos los instantes, hasta
𝐴𝑥𝑁−1 + 𝐵𝑢𝑁−1 − 𝑥𝑁 = 0 (48)
de modo que el conjunto resultante de ecuaciones sigue un patrón.
Ese sistema de ecuaciones puede escribirse de forma matricial como 𝐴𝑒𝑞𝑧 = 𝑏. Es exactamente la estructura de
las restricciones de igualdad. Por tanto, la matriz y el vector independiente se construyen como

𝐴𝑒𝑞 =
(

𝐵 −𝐼 0
0 𝐴 𝐵
0 0 0
0 0 0
−𝐼 0 0
𝐴 𝐵 −𝐼
⋱
0 0 0
−𝐼 0 0
𝐴 𝐵 −𝐼)

(49)

𝑏 =
(

−𝐴𝑥(𝑘)
0
0
⋮
0 )

(50)
Con ello, ya se encuentran redefinidos todos los objetos matemáticos que conforman el problema QP en base a
los valores, variables, parámetros y datos del MPC. El solver a elaborar será el encargado de resolver el problema
QP y estas matrices y vectores serán los datos o ingredientes que deben proporcionárseles para definir el
problema concreto que se esta resolviendo. Pues una vez definidos los ingredientes, es momento de ver las
instrucciones que el solver debe ejecutar para encontrar la solución óptima del problema.

3.2 Uso de ADMM para resolver MPC
Hasta aquí se han desarrollado las matemáticas necesarias para describir la técnica de control MPC con la
formulación típica de los problemas de optimización QP. Ahora, es necesario resolver ese problema QP y para
ello se ha contado con el algoritmo ADMM. Este es un algoritmo genérico para resolución de distintos tipos de
problemas de optimización. Como su alcance abarca más allá de problemas QP, la función objetivo general
suele venir descrita como (5)-(6).
Esta función (y sus restricciones) no se corresponde con la formulación usada para MPC una vez adaptada a QP.
Por tanto, es necesario hacer algunos ajustes más a la descripción del problema a la que se ha llegado en el
apartado anterior antes de aplicar ADMM. Así mismo, las instrucciones generales que conforman el algoritmo
(8)-(13) deben adaptarse a la formulación que está siendo usada.
3.2.1 Adaptación de la función objetivo
En primer lugar, resalta que la función objetivo propuesta (5) para ADMM está compuesta a su vez por dos
funciones distintas, cada una de ellas dependiente de una variable distinta. En consecuencia, la variable del
problema QP debe romperse en dos variables distintas y los elementos que componen la función objetivo deben
reagruparse en dos funciones individuales. Para la restricción propuesta en el enunciado del problema genérico
para ADMM (6) se seleccionan los valores de las matrices y el término independiente de forma que se pueda
reescribir la ecuación como

min
𝑧,𝑣
1
2
𝑧𝑇𝐻𝑧 + 𝑞𝑇𝑧 + 𝔗𝑒𝑞(𝑧) + 𝔗𝑖𝑛𝑒𝑞(𝑣) (51)
sujeto a
𝑧 = 𝑣 (52)
Así queda el problema original en dos subproblemas desacoplados a excepción de una ecuación lineal (52) que
los relaciona. Sin embargo, deja la incógnitade como expresar las restricciones de la formulación QP (35)-(36).
Estas restricciones se agregan como funciones indicadoras 𝔗𝑒𝑞(𝑧) y 𝔗𝑖𝑛𝑒𝑞(𝑣) a la función objetivo, del modo
que se expresa a continuación. La función objetivo queda entonces de la siguiente manera (51). Esta nueva
función objetivo está compuesta por dos funciones individuales en variables distintas (requerimiento de la
formulación para ADMM considerada).
Las funciones indicadoras siguen una determinada estructura. Se trata de funciones a trozos que solo toman dos
valores: cero e infinito. Valen cero cuando se cumplen la restricción que las generan (es decir, cuando la variable
pertenece al conjunto que define la ecuación de la restricción) e infinito en caso contrario. De este modo,
mientras que se cumplan las restricciones, las funciones indicadoras, que son un sumando más de la función
objetivo, no aportan nada. Esto es lo deseado, ya que mientras que se cumplan las restricciones, la función
objetivo se mantiene sin alterar. En caso de que no se cumpla la restricción, se está añadiendo un sumando de
valor infinito a la función objetivo, de modo que ese valor de las variables que incumple las restricciones jamás
será el valor que hace mínima la función. Las funciones indicadoras introducidas en el problema tienen la
estructura

𝔗𝑒𝑞 = {
0
+∞

si⁡𝐴𝑒𝑞𝑧 = 𝑏
en⁡otro⁡caso
(53)

𝔗𝑖𝑛𝑒𝑞 = {
0
+∞

si⁡𝑧 ≤ 𝑣 ≤ 𝑧
en⁡otro⁡caso
(54)
Volviendo a (51), la función de la variable z (tener en cuenta que ambas variables son iguales, es decir, la variable
de decisión), contiene tres términos. El primer de ellos es el término cuadrático de la formulación QP. El segundo
de ellos es el término lineal de la formulación original. El tercer término es la función indicadora que representa

Solver para MPC basado en ADMM

20
a la restricción de igualdad. Por otro lado, la función de la variable v, incluye solo a la función indicadora que
representa las restricciones de desigualdad.
En el apartado 2.3 se comentó como el método ADMM requiere el uso del lagraniano. El Lagrangiano resultante
de esta nueva reformulación del problema es

𝐿𝜌(𝑧, 𝑣, 𝜆) =
1
2
𝑧𝑇𝐻𝑧 + 𝑞𝑇𝑧 + 𝔗𝑒𝑞(𝑧) + 𝔗𝑖𝑛𝑒𝑞(𝑣) +
𝜌
2
‖𝑧 − 𝑣‖2 + 〈𝜆, 𝑧 − 𝑣〉 (55)
3.2.2 Desarrollo del algoritmo
Una vez el problema ha sido adaptado a la expresión deseada, es momento de aplicar el algoritmo ADMM. De
forma general, sigue la siguiente estructura (8)-(13). A continuación, se analiza cada una de las partes que lo
componen de forma individual. La condición inicial se mantiene igual, inicializando a cero los vectores
correspondientes. Lo siguiente son los tres pasos del algoritmo. El primero corresponde a la variable z y el
segundo a la variable v. En cada uno de ellos, se halla el valor de la variable que minimiza el Lagrangiano sujeto
a los valores de esa iteración para las otras. Es decir, la función a minimizar es dependiente solo de la variable
que se pretende calcular. Debido a eso, en cada paso se toma como función objetivo solo a la función individual
para esa variable más el par de sumandos que se añaden con el Lagrangiano, dando sentido a por qué separar la
función objetivo en dos funciones diferentes.
El resultado entonces para el primer paso del algoritmo es
1)
𝑧𝑘+1 = argmin
𝑧
1
2
𝑧𝑇𝐻𝑧 + 𝑞𝑇𝑧 +
𝜌
2
‖𝑧 − 𝑣𝑘‖
2
+ 〈(𝜆𝑘)
𝑇
, 𝑧〉
(56)
sujeto a
𝐴𝑒𝑞𝑧 = 𝑏
Puede observarse que la expresión vuelve a ser un problema de optimización en una única variable vectorial y,
por tanto, puede deshacerse la función indicadora para ser expresada como una restricción de igualdad de nuevo.
Si se desarrolla la norma y se agrupan términos, puede llegarse a la expresión
1
2
𝑧𝑇𝐻𝑧 + 𝑞𝑇𝑧 + (𝜆𝑘)
𝑇
𝑧 +
𝜌
2
𝑧𝑇𝐼𝑧 − 𝜌(𝑣𝑘)
𝑇
𝑧 +
𝜌
2
(𝑣𝑘)𝑇𝑣𝑘
=
1
2
𝑧𝑇(𝐻 + 𝜌𝐼)𝑧 + (𝑞 + 𝜆𝑘 − 𝜌𝑣𝑘)
𝑇
𝑧 =
1
2
𝑧𝑇𝐻𝑧𝑧 + (𝑞𝑧
𝑘)𝑇𝑧
(57)
Queda entonces el siguiente problema

𝑧𝑘+1 = argmin
𝑧
1
2
𝑧𝑇𝐻𝑧𝑧 + 〈(𝑞𝑧
𝑘)𝑇 , 𝑧〉
con⁡𝑞𝑧
𝑘 = (𝑞 + 𝜆𝑘 − 𝜌𝑣𝑘)
(58)
sujeto a
𝐴𝑒𝑞𝑧 = 𝑏
Intuitivamente se observa una estructura parecida con el problema QP originalmente planteado. La diferencia,
más allá del valor de los datos, es la ausencia de restricciones de desigualdad. Resulta que en esta situación sí
existe una solución analítica explicita al problema [18] dada por
𝑊𝜇 = −(𝐴𝑒𝑞𝐻𝑧
−1𝑞𝑧
𝑘 + 𝑏) con⁡𝑊 = (𝐴𝑒𝑞𝐻𝑧
−1𝐴𝑒𝑞
𝑇 ) (59)
𝑧𝑘+1 = −𝐻𝑧
−1(𝐴𝑒𝑞
𝑇 𝜇 + 𝑞𝑧
𝑘) (60)

Aplicando la solución, el cálculo de 𝑧𝑘+1 queda expresado como una única expresión que puede escribirse
fácilmente en código.
Para el segundo paso, puede aplicarse un desarrollo idéntico al primero, con lo que
2)
𝑣𝑘+1 = argmin
𝑣
𝜌
2
‖𝑧𝑘+1 − 𝑣‖
2
− 〈(𝜆𝑘)
𝑇
, 𝑣〉
(61) sujeto a
𝑧 ≤ 𝑣 ≤ 𝑧
pasa a ser

𝑣𝑘+1 = argmin
𝑧
1
2
𝑣𝑇(𝜌𝐼)𝑣 + 〈(𝑞𝑣
𝑘)𝑇 , 𝑣〉
con⁡𝑞𝑣
𝑘 = (−𝜆𝑘 − 𝜌𝑧𝑘+1)
(62)
sujeto a
𝑧 ≤ 𝑣 ≤ 𝑧
siguiendo
𝜌
2
𝑣𝑇𝐼𝑣 − 𝜌(𝑧𝑘+1)
𝑇
𝑣 +
𝜌
2
(𝑧𝑘+1)𝑇𝑧𝑘+1⁡−(𝜆𝑘)
𝑇
𝑣
=
1
2
𝑣𝑇(𝜌𝐼)𝑣 + (−𝜆𝑘 − 𝜌𝑧𝑘+1)
𝑇
𝑣 =
1
2
𝑣𝑇(𝜌𝐼)𝑣 + (𝑞𝑣
𝑘)𝑇𝑣
(63)
En este caso, la particularidad del problema resultante es que se encuentra totalmente desacoplado. Esto se debe
a la forma de la restricción y a que la norma no introduce acoplamientos. Debido al desacople, el problema
cuadrático en una variable vectorial (varias variables) puede descomponerse en un número de problemas
cuadráticos de una variable escalar (una única variable) independientes entre sí. El número de problemas
independientes es igual al tamaño de v. Los problemas cuadráticos en una única variable describen una parábola
en el plano. Es conocimiento matemático básico que el extremo relativo de una parábola se encuentra en el
vértice. Las restricciones de desigualdad representan un segmento en el que se observa la parábola, un intervalo
en el que debe encontrarse la solución. Por tanto, el mínimo absoluto debe ser el vértice o uno de los límites del
intervalo. El procedimiento queda entonces como sigue

𝑣𝑖
𝑘+1 = max [𝑧𝑖, min [𝑧𝑖,
−(𝑞𝑣
𝑘)
𝑖
𝜌
]] (64)
Se comprueba que valor es menor, si el vértice o la cota superior. Después, se toma el valor mayor entre la cota
inferior y el resultado de esta operación. Así, si el vértice se encuentra entre las dos cotas, el resultado es el
vértice. Si el vértice se encuentra por encima de la cota superior, el resultado será la cota superior, y si se
encuentra por debajo de la cota inferior, este será el resultado final.
El tercer y último paso de ADMM es (11), que tras aplicar la restricción planteada en (52), resulta en
3) 𝜆𝑘+1 = 𝜆𝑘 + 𝜌(𝑧𝑘+1 − 𝑣𝑘+1) (65)
Este paso no requiere mayor desarrollo, pues en su estado actual ya es implementable como código.
Finalmente, como fin del algoritmo, deben cumplirse unas condiciones de terminación o salida. La primera,
(12), asegura que el resultado óptimo encontrado cumple las restricciones propuestas. En el enunciado,
formulado de manera analítica, se esperaría que la expresión que representa las restricciones fuese igual a cero
como condición de salida. Sin embargo, al tratarse de un método numérico, es suficiente con que dicha expresión

Solver para MPC basado en ADMM

22
tome un valor menor a una determinada tolerancia, que puede ajustarse como parámetro de diseño del solver.
La segunda condición, (13), permite tener certeza de que se ha hallado un valor estable. Aplica el mismo
razonamiento sobre las tolerancias. Cabe añadir una tercera condición, de carácter computacional y no
relacionada con el problema matemático, que evite que el algoritmo se trabe en un bucle infinito. El resultado
de aplicar estas condiciones al problema exacto planteado es
‖𝑧𝑘+1 − 𝑣𝑘+1‖
2
≤ 𝜀𝑝 (66)
‖𝑧𝑘+1 − 𝑧𝑘‖
2
≤ 𝜀𝑑 (67)
𝑘 ≥ 𝑘𝑚𝑎𝑥 (68)
Una vez que se cumplenesas condiciones, el algoritmo ha finalizado. El valor de la variable 𝑧𝑘+1 es el valor
óptimo para el que la función objetivo es mínima. Por la construcción del vector de variables mostrada en el
apartado anterior, las primeras 𝑚 componentes, siendo 𝑚 el número de señales de control, son la acción de
control que se debe aplicar al sistema durante ese tiempo de muestreo para que se alcance la referencia.
3.3 Álgebra sparsa aplicada a ADMM
El desarrollo de solvers para problemas de optimización es un campo sobradamente estudiado. Como aporte de
este proyecto, se ha pretendido tomar esos conocimientos y aplicarlos al desarrollo en un lenguaje novedoso
como es Julia. En ese sentido se tiene que la eficiencia es un aspecto central de este proyecto, tal y como se ha
comentado en los apartados introductorios, y por ello no se ha limitado su búsqueda a cuestiones de lenguaje.
Desde el punto de vista del propio desarrollo, también pueden tomarse medidas para reducir el coste
computacional de todo el aparato matemático que se está manejando. Durante la aplicación del algoritmo
ADMM, se realizan numerosas operaciones con vectores y matrices. En base al número de estados y de señales
de control del sistema al que se enfrente el solver, estas entidades algebraicas pueden llegar a tener unas
dimensiones considerables. Es más, se da el caso de que algunas de estas matrices cuenten con gran cantidad de
componentes nulas. Así es el caso de la matriz 𝐻, que es diagonal a bloques. Estos espacios nulos ocupan gran
cantidad de memoria y suponen una serie de cálculos que no proporciona información útil sobre el problema.
Para ahorrar esfuerzo computacional en la ejecución del algoritmo ADMM se ha recurrido a la notación sparsa
en puntos claves del mismo.
El principal conjunto de operaciones matriciales que realizan en el desarrollo de ADMM es (59)-(60). El par de
ecuaciones puede ser expresado como
𝑊𝜇 = �̂�
(69)
con⁡�̂� = 𝐴𝑒𝑞𝐻𝑧
−1𝑞𝑧
𝑘 + 𝑏
−𝐻𝑧𝑧
𝑘+1 = �̂�
(70)
con⁡�̂� = 𝐴𝑒𝑞
𝑇 𝜇 + 𝑞𝑧
𝑘
Puede observarse que, en realidad, se trata de un par de sistemas de ecuaciones de la forma 𝐴𝑥 = 𝑏. Además,
los términos independientes están formados por una serie de operaciones que pueden reducirse a un producto
matriz-vector al que se le suma otro vector. La inversa de 𝐻𝑧, al ser esta diagonal a bloques y definida positiva,
sigue la misma estructura, por lo que convertirla a CSR es muy útil para resolver (70). La descomposición LDL
de 𝑊 también es muy sparsa y útil para resolver (69). Por ello, tanto la pareja de productos matriz-vector como
ambos sistemas de ecuaciones pueden resolverse de manera sparsa. Esto conlleva que los datos que se
proporcionen al solver, tanto como las operaciones realizadas, requieran modificarse de cara a implementar este
par de ecuaciones.

3.4 Implementación de un controlador MPC en bucle cerrado
De todo este conjunto de operaciones es de lo que se componen las entrañas del controlador. Pero hablar de un
controlador aislado tiene poco sentido. Es necesario colocarlo en su contexto y ver como se relaciona con los
demás participantes en el sistema. En especial, uno de los conceptos más importantes a tener en cuenta en la
estructura del sistema es la realimentación, independientemente de la técnica de control utilizada. La
realimentación es omnipresente en los sistemas de control, pues permite al sistema controlado adaptarse a
cambios en la referencia o rechazar los efectos de las perturbaciones. Para este proyecto, al estar centrado en la
técnica de control principalmente, se ha decidido no implantar una estructura de control muy complicada. Como
ya se ha mencionado, se ha descartado observar las salidas del sistema y lo que se pretende controlar son los
estados del sistema.
Se ha establecido que el sistema de control está compuesto simplemente por el controlador que implementa
MPC y por el sistema físico a controlar. El controlador recibe como entradas la referencia (estados y señales de
control) y calcula la señal de control adecuada, la óptima, obtenida como solución del problema QP. El
controlador entrega las señales de control al sistema, que modifica su estado en consecuencia. La única
interacción más que aparece en el sistema es la realimentación. Esta es representada por un flujo de información
del sistema al controlador. El controlador por MPC debe ser capaz de leer el estado actual del sistema en cada
tiempo de muestreo debido a que lo necesita para realizar las operaciones del cálculo de la señal de control
óptima.

Ilustración 2: Diagrama de control en bucle cerrado para MPC
Para poder testear el solver en funcionamiento, se debe realizar una simulación del sistema. Para ello, se enfrenta
al controlador a un modelo del sistema. Partiendo de unos valores iniciales, se simulan a través de un bucle los
ciclos de muestreo. En cada tiempo de muestro, actualizan los valores de 𝑥(𝑘) y (𝑥𝑟, 𝑢𝑟) en primer lugar. Luego
se resuelve el MPC. A continuación, se simula el sistema realizando las operaciones del modelo. No es
estrictamente necesario para el funcionamiento del sistema, pero, finalmente, puede guardarse los valores de las
variables que se consideren, para su posterior análisis. La Tabla 1 enseña el pseudocódigo que ejemplifica este
funcionamiento.
Tabla 1. Pseudocódigo de simulación
1 𝑥(𝑘) = 𝑥0
2 𝐟𝐨𝐫⁡1⁡𝐭𝐨⁡n_iteraciones
3 Actualizar⁡(𝑥𝑟 , 𝑢𝑟)
4 Resolver⁡MPC⁡ → 𝑢0
∗
5 Simular⁡sistema⁡𝑥(𝑘 + 1) = 𝐴𝑥(𝑘) + 𝐵𝑢0
∗
6 (Guardar⁡valores⁡𝑢(𝑘), 𝑧∗, 𝑘)
7 𝐞𝐧𝐝⁡𝐟𝐨𝐫

4 DESARROLLO DEL SOLVER

El solver constituye el elemento central del proyecto. Aunque el objetivo final se limite a crear el solver y
comprobar su correcto funcionamiento en simulación, el enfoque desde el que se realiza es que fuese ejecutable,
de manera eficiente, por un dispositivo electrónico capaz de controlar un sistema físico. Se hace especial énfasis
en la eficiencia para que la gama de dispositivos que pudiesen ejecutar el código abarque los dispositivos más
modestos posibles. Esta característica marca el diseño del mismo desde su propia estructura, que es presentada
a continuación.
En este apartado se realiza una descripción del conjunto de códigos que conforman la totalidad del solver, la
estructura que siguen y su diseño. Se explican las decisiones de diseño tomadas y las dificultades afrontadas
durante el desarrollo. De las explicaciones se deducen y se exponen los pseudocódigos de las funciones
implementadas. El código final implementado en los respectivos lenguajes se adjunta en los anexos de este
documento.
4.1 Desarrollo en Matlab
Una de las decisiones de diseño más importantes es el lenguaje en el que debe ser implementado. Como ya se
ha explicado, la elección tomada fue Julia. Sin embargo, debido a la cantidad de herramientas que ofrece y la
potencia de su entorno de desarrollo, se decidió hacer una versión preliminar en Matlab. Esto ha permitido
familiarizarse de forma más cómoda con la idiosincrasia de MPC. Una vez escrito todo el código en Matlab, ha
sido transcrito a Julia haciendo los ajustes correspondientes.
4.1.1 Descripción del solver
Un solver es un programa que aplica una serie de operaciones matemáticas sobre un conjunto de datos a través
de un algoritmo para proporcionar la solución a un determinado problema matemático caracterizado por ese
conjunto datos. En el contexto de este proyecto ese problema matemático es un problema de optimización tipo
QP. Como el problema de optimización concreto deriva de la descripción del sistema de control MPC,
recontextualizado, implica que es necesario distinguir entre dos bloques diferenciados entre sí para la
implementación: la adecuación del problema de control al problema matemático y la solución del problema
matemático en sí.
Esta es la principal diferenciación