Rápido
En este artículo lo haremos:
1. Introducción
El control predictivo del modelo (MPC) es una metodología de control de comentarios popular en la que un problema de control óptimo de horizonte finito (OCP) se resuelve iterativamente con un estado medido actualizado en cada iteración.
En el OCP se usa un modelo de la planta para encontrar el control óptimo de bucle abierto sobre el horizonte finito considerado. Debido a que el modelo no puede capturar el comportamiento de la planta verdadera al 100%, y porque en el mundo real el sistema está sujeto a ruido y perturbaciones, uno solo aplica la primera parte del control óptimo de bucle abierto, y regularmente medidas el estado para volver a solucionar el OCP. Esto cierra el bucle y crea un comentario.
Las matemáticas detrás de esto son relativamente simples e intuitivas (especialmente en comparación con cosas como el control robusto) y es fácil codificar un controlador MPC. Otros profesionales son que puede manejar las limitaciones duras y suaves en el estado y el control (las limitaciones duras son estrictas, mientras que las restricciones suaves se aplican mediante la penalización en la función de costo) y generalmente se pueden usar en sistemas no lineales con restricciones no convexas (dependiendo de cuán desagradables son, por supuesto!). La única estafa es que necesita resolver los problemas de optimización “en línea” en tiempo real, lo que puede ser un problema si está controlando un sistema rápido o tiene recursos informáticos limitados.
1.2 Ejemplo de ejecución
A lo largo del artículo, consideraré el doble integrador con un control de retención de orden cero (es decir, un control constante por partes) como ejemplo en ejecución en el código. El sistema de tiempo continuo dice:
\[
\begin{align*}
\dot x _1
con \ (t \ in \ mathbb {r} \) denotando tiempo. Aquí \ (x_1
Si limitamos que el control sea constante por partes en intervalos de longitud 0.1 segundos, obtenemos el sistema de tiempo discreto:
\[
\begin{align*}
x_{k+1} &= A x_k + Bu_k,
\end{align*}
\]
con \ (k \ in \ mathbb {z} \), donde,
\[
\begin{align*}
A :=
\left(
\begin{array}{cc}
1 & 0.1\\
0 & 1
\end{array}
\right), \,\,
B :=
\left(
\begin{array}{c}
0.05\\
1
\end{array}
\right)
\end{align*}
\]
y \ (x_k \ in \ mathbb {r}^2 \), \ (u_k \ in \ mathbb {r} \). (Tenga en cuenta que uso \ (x_k \) y \ (u_k \) para referirme al estado y control del sistema de tiempo discreto para hacer que la notación sea más nater en el resto del artículo).
Puede usar la función Cont2DIScrete de los paquetes Scipy para obtener este sistema de tiempo discreto, de la siguiente manera:
import numpy as np
from scipy.signal import cont2discrete
A = np.array([[0, 1],[0, 0]])
B = np.array([[0],[1]])
C = np.array([[1, 0],[0, 1]])
D = np.array([[0, 0],[0, 0]])
dt = 0.1 # in seconds
discrete_system = cont2discrete((A, B, C, D), dt, method='zoh')
A_discrete, B_discrete, *_ = discrete_system2. Problema de control óptimo
Consideraremos el siguiente problema de control óptimo de tiempo discreto (OCP):
\[
\begin{equation}
{\mathrm{OCP}(\bar x):}\begin{cases}
\min\limits_{\mathbf{u},\mathbf{x}} \quad & \sum_{k=0}^{K-1} \left(x_k^{\top}Qx_k + u_k^\top R u_k \right)+ x_{K}^\top Q_K x_{K}\\
\mathrm{subject\ to:}\quad & x_{k+1} = Ax_k + Bu_k, &k \in[0:K-1]& \ dots (1) \\
\ quad & x_0 = \ bar x, & & \ dots (2) \\
\ quad & x_k \ in [-1,1]\ Times (-\ Infty, \ Infty), & k \ in[1:K]& \ dots (3) \\
\ quad & u_k \ in [-1,1],& familiares[0:K-1]& \ dots (4)
\ end {casos}
\ end {ecuación}
\]
dónde,
- \ (K \ in \ mathbb {r} _ {\ geq 0} \) denota el horizonte finito sobre el cual resolvemos el OCP,
- \ (k \ in \ mathbb {z} \) denota el paso de tiempo discreto,
- \ ([p:q]\), con \ (p, q \ in \ mathbb {z} \), denota el conjunto de enteros \ (\ {p, p+1, …, q \} \),,
- \ (\ bar x \ in \ mathbb {r}^2 \) denota el condición inicial del sistema dinámico,
- \ (x_k \ in \ mathbb {r}^2 \) denota el estado en el paso \ (k \),
- \ (u \ in \ mathbb {r} \) denota el control en el paso \ (k \),
- \ (Q \ in \ mathbb {r}^{2 \ times 2} \), \ (r \ in \ mathbb {r} \) y \ (q_k \ in \ mathbb {r}^{2 \ times 2} \) son matrices definidas cuadradas que especifican la función de costo (\ (r \) es escalar aquí porque \ (u \) es escalar).
Además, lo dejaremos
- \ (\ Mathbf {u}: = (u_0, u_1,…, u_ {k – 1}) \ in \ mathbb {r}^k \) denota el secuencia de control,
- \ (\ \ mathbf {x}: = (x_0, x_1, …, x_ {k}) \ in \ mathbb {r}^{2 (k+1)} \) denota el Secuencia de estado.
Para ser riguroso, diremos que el par \ ((\ mathbf {u}^{*}, \ mathbf {x}^{*}) \ in \ mathbb {r}^k \ times \ mathbb {r}^{2 (k+1)} \) es un solución a \ (\ mathrm {OCP} (\ bar {x}) \) siempre que minimice la función de costo en todos los pares admisibles, es decir, es,
\[
\begin{equation*}
J(\mathbf{u}^{*}, \mathbf{x}^{*}) \leq J(\mathbf{u}, \mathbf{x}),\quad \forall (\mathbf{u},\mathbf{x})\in \Omega
\end{equation*}
\]
donde \ (j: \ mathbb {r}^k \ times \ mathbb {r}^{2 (k+1)} \ rectarrow \ mathbb {r} _ {\ geq 0} \) es, es,
\[
\begin{equation*}
J(\mathbf{u},\mathbf{x}) :=\left( \sum_{k=0}^{K-1 }x_k^\top Q x_k + u_k^\top R u_k \right) + x_K^\top Q_K x_K
\end{equation*}
\]
y \ (\ omega \) denota todos los pares admisibles,
\[
\Omega :=\{(\mathbf{u},\mathbf{x})\in \mathbb{R}^{K}\times \mathbb{R}^{2(K+1)} : (1)-(4)\,\, \mathrm{hold}\}.
\]
Por lo tanto, el problema de control óptimo es encontrar una secuencia de control y estado, \ ((\ mathbf {u}^{*}, \ mathbf {x}^{*}) \), que minimiza la función de costo, sujeto a la dinámica, \ (f \), así como restricciones en el estado y el control, \ (x_k \ in in in \ in in in in in in \ in in in \ in in in \ in in \ in in \ in in \ in in in \ in in \ in in in \ in in \ in in in \ in in \ in in in \ in in in \ in in \ in in \ in in in \ in in \ in in in \ in in \ in in in \ in in \ in in \ in in \ in \ in in in \ in in \ in \ in \ in in in \ in in en in[-1,1]\ Times (-\ Infty, \ Infty) \), \ (u_k \ in [-1,1] \), para todos \ (k \). La función de costo es vital para el rendimiento del controlador. No solo en el sentido de garantizar que el controlador se comporte bien (por ejemplo, para evitar señales erráticas) sino también para especificar el punto de equilibrio El estado de circuito cerrado se asentará. Más sobre esto en la Sección 4.
Tenga en cuenta cómo \ (\ mathrm {OCP} (\ bar x) \) se parametriza con respecto al estado inicial, \ (\ bar x \). Esto proviene de la idea fundamental detrás de MPC: que el problema de control óptimo se resuelve iterativamente con un estado medido actualizado.
2.1 Codificación de un solucionador OCP
Casadi’s opti Stack hace que sea realmente fácil configurar y resolver el OCP.
Primero, algunos preliminares:
from casadi import *
n = 2 # state dimension
m = 1 # control dimension
K = 100 # prediction horizon
# an arbitrary initial state
x_bar = np.array([[0.5],[0.5]]) # 2 x 1 vector
# Linear cost matrices (we'll just use identities)
Q = np.array([[1. , 0],
[0. , 1. ]])
R = np.array([[1]])
Q_K = Q
# Constraints for all k
u_max = 1
x_1_max = 1
x_1_min = -1Ahora definimos las variables de decisión del problema:
opti = Opti()
x_tot = opti.variable(n, K+1) # State trajectory
u_tot = opti.variable(m, K) # Control trajectoryA continuación, imponemos las restricciones dinámicas y configuramos la función de costo:
# Specify the initial condition
opti.subject_to(x_tot[:, 0] == x_bar)
cost = 0
for k in range(K):
# add dynamic constraints
x_tot_next = get_x_next_linear(x_tot[:, k], u_tot[:, k])
opti.subject_to(x_tot[:, k+1] == x_tot_next)
# add to the cost
cost += mtimes([x_tot[:,k].T, Q, x_tot[:,k]]) + \
mtimes([u_tot[:,k].T, R, u_tot[:,k]])
# terminal cost
cost += mtimes([x_tot[:,K].T, Q_K, x_tot[:,K]])def get_x_next_linear(x, u):
# Linear system
A = np.array([[1. , 0.1],
[0. , 1. ]])
B = np.array([[0.005],
[0.1 ]])
return mtimes(A, x) + mtimes(B, u)El código mTimes ([x_tot[:,k].T, q, x_tot[:,k]]) implementa la multiplicación de matriz, \ (x_k^{\ top} q x_k \).
Ahora agregamos las limitaciones de control y estado,
# constrain the control
opti.subject_to(opti.bounded(-u_max, u_tot, u_max))
# constrain the position only
opti.subject_to(opti.bounded(x_1_min, x_tot[0,:], x_1_max))y resolver:
# Say we want to minimise the cost and specify the solver (ipopt)
opts = {"ipopt.print_level": 0, "print_time": 0}
opti.minimize(cost)
opti.solver("ipopt", opts)
solution = opti.solve()
# Get solution
x_opt = solution.value(x_tot)
u_opt = solution.value(u_tot)Podemos trazar la solución con la función trapt_solution () del repositorio.
from MPC_tutorial import plot_solution
plot_solution(x_opt, u_opt.reshape(1,-1)) # must reshape u_opt to (1,K)3. Control predictivo del modelo
La solución a \ (\ mathrm {OCP} (\ bar x) \), \ ((\ mathbf {x}^{*}, \ mathbf {u}^{*}) \), para un dado \ (\ bar x \), nos proporciona un bucle abierto Control, \ (\ mathbf {u}^{*} \). Nosotros ahora Cerrar el bucle resolviendo iterativamente \ (\ Mathrm {OCP} (\ Bar X) \) y actualizando el estado inicial (este es el algoritmo MPC).
\[
\begin{aligned}
&\textbf{Input:} \quad \mathbf{x}^{\mathrm{init}}\in\mathbb{R}^2 \\
&\quad \bar x \gets \mathbf{x}^{\mathrm{init}} \\
&\textbf{for } k \in [0:\infty) \textbf{:} \\
&\quad (\mathbf{x}^{*}, \mathbf{u}^{*}) \gets \arg\min \mathrm{OCP}(\bar x)\\
&\quad \mathrm{apply}\ u_0^{*} \mathrm{\ to\ the\ system} \\
&\quad \bar x \gets \mathrm{measured \ state\ at\ } k+1 \\
&\textbf{end for}
\end{aligned}
\]
3.1 Codificación del algoritmo MPC
El resto es bastante sencillo. Primero, pondré todo nuestro código anterior en una función:
def solve_OCP(x_bar, K):
....
return x_opt, u_optTenga en cuenta que he parametrizado la función con el estado inicial, \ (\ bar x \), así como el horizonte de predicción, \ (k \). El bucle MPC se implementa luego con:
x_init = np.array([[0.5],[0.5]]) # 2 x 1 vector
K = 10
number_of_iterations = 150 # must of course be finite!
# matrices of zeros with the correct sizes to store the closed loop
u_cl = np.zeros((1, number_of_iterations))
x_cl = np.zeros((2, number_of_iterations + 1))
x_cl[:, 0] = x_init[:, 0]
x_bar = x_init
for i in range(number_of_iterations):
_, u_opt = solve_OCP(x_bar, K)
u_opt_first_element = u_opt[0]
# save closed loop x and u
u_cl[:, i] = u_opt_first_element
x_cl[:, i+1] = np.squeeze(get_x_next_linear(x_bar,
u_opt_first_element))
# update initial state
x_bar = get_x_next_linear(x_bar, u_opt_first_element)Nuevamente, podemos trazar la solución de circuito cerrado.
plot_solution(x_cl, u_cl)Tenga en cuenta que he “medido” el estado de la planta utilizando la función get_x_next_linear (). En otras palabras, he asumido que nuestro modelo es 100% correcto.
Aquí hay una parcela del circuito cerrado de un montón de estados iniciales.
4. Otros temas
4.1 Estabilidad y viabilidad recursiva
Dos de los aspectos más importantes de un controlador MPC son viabilidad recursiva del OCP iterativamente invocado y estabilidad del circuito cerrado. En otras palabras, si he resuelto el OCP en el momento \ (k \), ¿existirá una solución al OCP en el momento \ (k+1 \)? Si existe una solución al OCP en cada paso de tiempo, ¿el estado de circuito cerrado se asentará asintóticamente en un punto de equilibrio (es decir, será estable)?
Asegurar que un controlador MPC exhiba estas dos propiedades implica el diseño cuidadosamente de la función de costo y las limitaciones, y la elección de un horizonte de predicción lo suficientemente largo. Volviendo a nuestro ejemplo, recuerde que las matrices en la función de costo simplemente fueron elegidas para ser:
\[
Q = \left( \begin{array}{cc}
1 & 0\\
0 & 1
\end{array}
\right),\, Q_K = \left( \begin{array}{cc}
1 & 0\\
0 & 1
\end{array}
\right),\, R = 1.
\]
En otras palabras, el OCP penaliza la distancia del estado al origen y, por lo tanto, lo impulsa allí. Como probablemente pueda adivinar, si el horizonte de predicción, \ (k \), es muy corto y el estado inicial se encuentra muy cerca de las restricciones en \ (x_1 = \ pm 1 \), el OCP encontrará soluciones con “previsión” insuficiente y el problema será infalible en alguna iteración futura del bucle MPC. (También puede experimentar con esto haciendo \ (k \) pequeño en el código).
4.2 Algunos otros temas
MPC es un campo de investigación activo y hay muchos temas interesantes para explorar.
¿Qué pasa si no se puede medir el estado completo? Esto se relaciona con observabilidad y Salida MPC. ¿Qué pasa si no me importa la estabilidad asintótica? Esto (a menudo) tiene que ver con MPC económico. ¿Cómo hago el controlador? robusto al ruido y perturbaciones? Hay algunas maneras de lidiar con esto, con Tubo MPC Probablemente siendo el más conocido.
Los artículos futuros podrían centrarse en algunos de estos temas.
5. Lectura adicional
Aquí hay algunos libros de texto estándar y muy buenos en MPC.
[1] Grüne, L. y Pannek, J. (2016). Control predictivo del modelo no lineal.
[2] Rawlings, JB, Mayne, DQ y Diehl, M. (2020). Control predictivo del modelo: teoría, cálculo y diseño.
[3] Kouvaritakis, B. y Cannon, M. (2016). Control predictivo del modelo: clásico, robusto y estocástico.