Apuntes de Control Predictivo

William Colmenares
5 de febrero de 2003

Índice
1. Predicciones en Variables de Estado 1

2. Formulación Predictiva 3

3. Una formulación alternativa 5

4. Obervadores de estado 7

5. Manejo de incertidumbre y perturbaciones 8

1. Predicciones en Variables de Estado

Recordemos que cuando hablamos de Control Predictivo lo fundamental es
determinar las predicciones en función del modelo del que disponemos.
En esta sección nos ocuparemos de la formulación explı́cita de las predic-
ciones para el caso en el que tengamos una representación en Variables de Estado
del sistema a controlar. Esto es:
x(t + 1) = Ax(t) + Bu(t)
ym (t) = Cx(t); x(0) = x0 , (1)
donde x ∈ Rn es el vector de estados y u ∈ Rm es el vector de control y
y ∈ Rp es el vector de salida . A, B, C son matrices constantes de dimensiones
adecuadas.
A partir de 1 es fácil formular las predicciones hasta el horizonte de predic-
ción, N . Para simplificar la descripción, asumiremos que t = 0, esto es, en el
instante actual calcularemos la ley de control.

x(1) = Ax(0) + Bu(0)

x(2) = A2 x(0) + ABu(0) + Bu(1)
x(3) = A3 x(0) + A2 Bu(0) + ABu(1) + Bu(2) (2)
..
.
x(N ) = AN x(0) + AN −1 Bu(0) + . . . + ABu(N − 2) + Bu(N − 1)

1
 Observe que para poder hacer las predicciones sólo se requiere la condición
inicial x(0) y el vector de controles futuros, u(0), . . . , u(N −1). El Control Predic-
tivo se ocupa de calcular esos controles futuros que minimizan alguna función de
costo, que precisaremos más adelante, aplicando solamente, el primero de ellos,
u(0).
En forma matricial, 2, resulta:
      
x(1) B 0 ... 0 u(0) A
 x(2)   AB B 0 ... 0   u(1)   A2 
      3 
 x(3)   A2 B AB B 0 ... 0    u(2)   A 
 = +  x(0)
 ..   ..   ..   .. 
 .   .   .   . 
x(N ) AN −1 B AN −2 B ... B u(N − 1) AN
(3)
Normalmente, las predicciones se realizan, no con el vector de control u(k),
k = 0, . . . , N − 1, sino con su incremento, ∆u(k) = u(k) − u(k − 1). Ello para
introducir un integrador en la estrategia de Control Predictivo, con los beneficios
de eliminación de error estacionario o rechazo a perturbaciones del mismo tipo
en algunos sistemas, que ello supone. Observemos que:

u(0) = ∆u(0) + u(−1)

u(1) = ∆u(1) + ∆u(0) + u(−1)
..
.
u(N − 1) = ∆u(N − 1) + ∆u(N − 2) + . . . + ∆u(0) + u(−1), (4)

donde u(−1) es la última ley de control aplicada.

Al reemplazar en 3 el incremento, resulta:
    
x(1) B 0 ... 0 ∆u(0)
 x(2)   AB + B B 0 ... 0  ∆u(1) 
    
 x(3)   A2 B + AB + B AB + B B ... 0  ∆u(2) 
 =  +
 ..   ..  .. 
 .   .  . 
PN N −j
PN N −j
x(N ) A
j=1  B j=2 A  B . . .  B ∆u(N − 1)

B A
 AB + B   A2 
 2   3 
   
+  A B + AB + B  u(−1) +  A  x(0)
 ..   .. 
 .   . 
PN
j=1 A
N −j
B AN
(5)
Observe que para poder hacer las predicciones de los estados, x(k), k =
1, . . . , N , adicional a la condición inicial x(0), se requiere de la última ley de
control aplicada u(−1).
Comentamos que en muchas circunstancias, como por ejemplo sistemas con
respuesta inversa o con retardo, sólo se requieren las predicciones a partir de un
cierto instante T > 0. De la misma forma, para simplificar el problema y hacerlo

2
menos pesado en el cálculo de las estrategias óptimas de control, el número de
cambios posibles se reduce de N a Nu < N (horizonte de control más pequeño
que el de predicción). Ello se traduce en la eliminación de algunas filas (primer
caso) o columnas (segundo caso) de las matrices de datos en 5.

2. Formulación Predictiva
Con las predicciones de la sección anterior, ya podemos abordar el problema
de control óptimo que, en cada iteración, resuelve el enfoque predictivo.
De esta forma, para un sistema como el descrito por 1, en cada iteración,
y desde una condición inicial x(0) dada, en el enfoque predictivo se calcula
la secuencia de controladores que se aplicarı́an al sistema, que minimizan la
siguiente función de costo:
N
X Nu
¡ ¢T ¡ ¢ X
J(y, u) = y(j) − yref (j) Qj y(j) − yref (j) + ∆u(j)T Rj ∆u(j) (6)
j=1 j=0

En 6, Qj y Rj son matrices conocidas semidefinidas positivas (esto implica

que la función J(y, u) ≥ 0). Observe que el horizonte de predicción N no tiene
que ser el mismo que el de control Nu y que por simplicidad, hemos fijado al
horizonte de predicción comenzando en el instante siguiente al presente. yref es
la trayectoria de referencia que deseamos que el sistema “siga”.
Se desea minimizar 6 sujeto a 1 y normalmente se imponen condiciones
adicionales del tipo:
ymin ≤ y(k) ≤ ymax
umin ≤ u(q) ≤ umin (7)
∆umin ≤ ∆u(q) ≤ ∆umax
para todo k = 1, . . . , N y q = 0, . . . , Nu .
Antes de calcular el controlador predictivo, hagamos las siguientes defini-
ciones:    
x(1) ∆u(0)
 x(2)   ∆u(1) 
   
 x(3)   ∆u(2) 
X=  ; ∆U =  .
 ..   .. 
 .   . 
x(N ) ∆u(N − 1)
Nótese que ∆U es el vector de leyes de control futuras que queremos deter-
minar (variables de decisión).
Sean igualmente:  
y(1)
 y(2) 
 
 
Y =  y(3)  = CX, (8)
 .. 
 . 
y(N )

3
donde  
C 0 ... 0
 0 C 0 ... 0 
 
 0 0 C 0 ... 0 
C= . (9)
 .. 
 . 
0 ... 0 C
Podemos ahora escribir las ecuaciones de predicción de la salida en forma
compacta (matricial), utilizando 5 y 8, resultando:
¡ ¢
Y = C Ax(0) + Bu(−1) + Θ∆U . (10)
Es fácil extraer de 5 la definición de las matrices A, B y Θ en 10.
A partir de 10 podemos escribir en forma matricial el problema de opti-
mización que resolvemos en control predictivo en cada iteración. A saber, el
funcional 6 se describe como:
J(∆U ) = (Y − Yref )T Q(Y − Yref ) + ∆UT R∆U (11)
donde:
     
yref (1) Q1 0 ... 0 R1 0 ... 0
 yref (2)   0 Q2 0 ... 0   0 R2 0 ... 0 
     
 yref (3)   0 0 Q3 0 ... 0   0 0 R3 0 ... 0 
Yref = ;Q =  ;R =  
 ..   ..   .. 
 .   .   . 
yref (N ) 0 ... 0 QN 0 ... 0 RN
Sin pérdida de generalidad, en lo que sigue supondremos que Q = I (matriz
identidad) y R = λ (escalar).
El problema de Control Predictivo se reduce a:
¡ ¢ ¡ ¢T
mı́n ∆UT ΘT C T CΘ + λI ∆U + 2 CAx(0) + CBu(−1) − Yref Θ∆U (12)
∆U

sujeto a:
restricciones en las salidas
Ymin − CAx(0) − CBu(−1) ≤ CΘ∆U ≤ Ymax − CAx(0) − CBu(−1)

restricciones en los saltos en el control

∆Umin ≤ ∆U ≤ ∆Umax

restricciones en el valor absoluto del control

 
Im 0 . . . 0
 Im Im 0 . . . 0 
 
 Im Im Im 0 . . . 0 
Umin − U−1 ≤   ∆U ≤ Umax − U−1
 .. 
 . 
Im Im Im Im . . . Im

4
 Ymin , Ymax , U−1 , Umax ,Umin ,∆Umin ,∆Umax son los vectores que con-
tienen (en stack), los valores lı́mite y el último control aplicado. Im es la matriz
identidad de dimensión m × m.

3. Una formulación alternativa

El usar 5 como predictor presenta el inconveniente de que para horizontes de
predicción grandes, como es el caso usual en cualquier estrategia predictiva en el
que N = 20 ó más, se pueden introducir problemas numéricos cuando la matriz
del sistema (A) no esté bien condicionada, i.e., cuando haya mucha diferencia
entre su autovalor más grande y el más pequeño. Ello, porque para predecir hay
que elevar esa matriz a potencias del orden del horizonte de predicción. Para
solventar ese problema, se incluyen como variables de decisión (incógnitas), la
predicción de los estados del sistema. Ello simplifica bastante la formulación del
problema y elimina las potencias de la matriz A, mejorando de esta manera la
robustez numérica del algoritmo; se aumentan considerablemente, sin embargo,
el número de variables de decisión y se introducen restricciones de igualdad que
antes no aparecen en la formulación original.
Para esta nueva formulación recordemos que:

x(1) = Ax(0) + Bu(0)

x(2) = Ax(1) + Bu(1)
x(3) = Ax(2) + Bu(2)
..
. (13)
x(N ) = Ax(N − 1) + Bu(N − 1)

Si sustituimos el valor absoluto u por su valor relativo ∆u (ver 4), obtenemos:

x(1) − B∆u(0) = Ax(0) + Bu(−1)

¡ ¢
x(2) − Ax(1) − B ∆u(1) + ∆u(0) = Bu(−1)
¡ ¢
x(3) − Ax(2) − B ∆u(2) + ∆u(1) + ∆u(0) = Bu(−1)
..
. (14)
¡ ¢
x(N ) − Ax(N − 1) − B ∆u(N − 1) + . . . + ∆u(0) = Bu(−1)

Antes de escribir la formulación matricial (compacta) de la estrategia alter-

nativa, permı́tasenos definir las siguientes matrices:
   
I 0 ... 0 B 0 ... 0
 −A I 0 . . . 0   B B 0 ... 0 
   
 ... 0   
α= 0 −A I 0  ;β =  B B B 0 ... 0 
 ..   .. 
 .   . 
0 ... 0 −A I (N ×n)×(N ×n)
B B ... B (N ×n)×(N ×m)

5
 Igualmente, definamos en nuevo vector de variables de decisión,
 
x(1)
 x(2) 
 
 x(3) 
 
 .. 
 . 
 
 x(N ) 

Vd =  .
∆u(0) 
 
 ∆u(1) 
 
 ∆u(2) 
 
 .. 
 . 
∆u(N − 1)
y sean (con relación a 9):
£ ¤ £ ¤
C1 = C 0 ¡ ¢ ¡ ¢ ; C2 = 0 IN ×m ¡ ¢ ¡ ¢,
N ×p × N ×(n+m) N ×p × N ×(n+m)

de modo que Y = C1 Vd y ∆U = C2 Vd .
Las ecuaciones para predicción son:
   
A B
 0   B 
£ ¤    
   
α β Vd =  0  x(0) +  B  u(−1). (15)
 ..   .. 
 .   . 
0 B
Observe que en las ecuaciones de predicción no hay ahora potencias de la
matriz A.
El problema de control predictivo puede reescribirse, en el marco de la nueva
formulación de la forma:
T
mı́n VdT (C1T C1 + λC2T C2 )Vd − 2Yref C1 Vd (16)
Vd

sujeto a:
restricciones en la salida del sistema
Ymin ≤ C1 Vd ≤ Ymax

restricciones en la tasa de cambio del control

∆Umin ≤ C2 Vd ≤ ∆Umax

restricciones en el valor absoluto del control

 
Im 0 . . . 0
 Im Im 0 . . . 0 
 
 Im Im Im 0 . . . 0 
Umin − U−1 ≤   C2 Vd ≤ Umax − U−1
 .. 
 . 
Im Im Im Im . . . Im

6
 restricciones en las condiciones dinámicas expresadas en las igualdades de
15

4. Obervadores de estado
Cualquiera que sea la formulación del controlador predictivo, en ambas se
requiere de la condición inicial x(0). En muchos casos, no todos los estados
están disponibles en la salida del medible del sistema –no hay medición directa
en la salida de todos los estados–, y aunque los estados estuviesen disponibles,
los modelos con los que se hace la predicción son sólo aproximaciones del sis-
tema real y además se ven afectados por perturbaciones externas que no siempre
pueden medirse. En cualquiera de los casos anteriores, sólo se cuenta con infor-
mación parcial de los estados reales del sistema (aquella que puede inferirse en
la salida) y debe recurrirse a observadores para poder estimar su evolución.
Para desarrollar el observador, adoptemos la siguiente nomenclatura: ξ(k +
1|k) es el estimado de los estados en k + 1 dado que conocemos la salida en k.
ξ(k|k) es entonces el estimado de los estados en k dado que se conoce la salida
en k.
El observador, propuesto originalmente por Kalman, es el siguiente (con
relación a 1) :
¡ ¢
ξ(k|k) = ξ(k|k − 1) + L y(k) − Cξ(k|k − 1) (17)
ξ(k + 1|k) = Aξ(k|k) + Bu(k) (18)

Combinando las dos ecuaciones anteriores se obtiene:

ξ(k + 1|k) = (A − KC)ξ(k|k − 1) + Bu(k) + Ky(k) (19)

donde K = AL.
Si los autovalores de A − KC están en el cı́rculo unitario, el observador
converge asintóticamene al valor real del estado que queremos estimar.
Aunque hay sólidas estrategias para el cálculo de K ó L, acá sólo nos pre-
ocuparemos porque el observador ξ tenga una dinámica más rápida que la del
sistema original 1, tratando a su vez de que la ganancia K no sea muy grande
porque ello podrı́a introducir grandes oscilaciones en el estimador e inclusive
inestabilidad. Para estimar la ganancia entonces nos limitaremos a cualquier
método de ubicación de polos.
Observe que el observador es un sistema del mismo orden que el sistema
original y su salida es la condición inicial x(0) que usaremos en la próxima
iteración, si usamos 19, o la condición inicial x(0) que usaremos para calcular
u(0) en la iteración actual, si usamos 17. No hay mucha diferencia entre una y
otra y si se usa por ejemplo, Simulink, sólo se usarı́a la segunda por permitir
escribir una función de transferencia.

7
5. Manejo de incertidumbre y perturbaciones
Además de las ventajas numéricas que presenta la formulación en Variables
de Estado, el tratamiento de las incertidumbre en el modelo y las perturbaciones
externas no medibles (desde el punto de vista entrada salida, son indistingibles
la una de la otra en los sistemas lineales), es bastante más sencillo y natural
en la formulación de estados con observador. Otras estrategias, como DMC y
GPC, pueden considerarse como casos particulares del primer enfoque.
En lo que sigue cuando hacemos referencia a perturbación, nos referimos de
forma general a la incertidumbre en el modelo, a las perturbaciones externas
no medibles o a ambas. De hecho, consideraremos que afecta al sistema como
una perturbación externa que colocaremos comúnmente añadida a la salida del
modelo, ello sin pérdida de generalidad.
Para considerar la incertidumbre lo primero que debe considerarse es su na-
turaleza, determinı́stica o estocástica. En el primer caso, comunmente hacemos
un estudio comparativo entre la salida del sistema y del modelo y tratamos de
ajustar una curva conocida a la diferencia (sacar una tendencia) o, en muchos ca-
sos conocemos la naturaleza de los cambios que afectan al proceso (por ejemplo,
cambios en la carga), pero no se conocen a priori cuándo y con qué intensidad
ocurriran. En ete caso –determinı́stico– se conoce la forma de la señal pero no se
conocen sus parámetros. Ası́ puede ser un escalón (o salto) del que no se conoce
la altura, o una rampa de la que no se conoce pendiente y punto de inicio o una
sinusoide de la que se conoce frecuencia pero no amplitud y desfasaje. En este
caso sólo se requiere incorporar al sistema el generador de la señal conocida.
En el caso estocástico, por ejemplo el vaiven de las olas que golpean el
casco de tanquero, además de lo anterior, la señal es generada por una variable
estocástica cuyas propiedades se conocen.
La única diferencia entre ambos es que uno no posee señal de entrada que
genera la incertidumbre y el otro si. Para efectos del cálculo del control predictivo
da igual cualquier enfoque, no ası́ para calcular la salida del modelo que debe
incluir la señal estocástica externa.
Consideremos que la perturbación es determinı́stica. La salida del sistema
será:
y(k) = ym (k) + d(k) (20)
Observe que de no haber perturbación la salida del modelo y del sistema
coinciden.
Supongamos ahora que conocemos la expresión para la evolución de la per-
turbación, esto es:
ζ(k + 1) = Aζ ζ(k) (21)
d(k) = Cζ ζ(k) (22)
Considerar la perturbación equivale, simplemente, a escribir en representación
de estados, la ecuación 20, esto es:
x(k + 1) = Ax(k) + Bu(k)

8
 ζ(k + 1) = Aζ ζ(k) (23)
· ¸
£ ¤ x(k)
y(k) = C C ζ
ζ(k)

Recordemos ahora, las señales más comunes usadas en control predictivo.

Ellas son:
Perturbaciones constantes (escalón)

d(k + 1) = d(k)

Esto es Aζ = 1, Cζ = 1.
Rampa

x1 (k + 1) = x(2)
x2 (k + 1) = 2x(2) − x(1)
d(k) = x2 (k)

siendo · ¸
0 1
Aζ =
−1 2

Sinusoide (variaciones periódicas) con frecuencia ω0 y tiempo de muestreo

Ts conocidos

x1 (k + 1) = x(2)
x2 (k + 1) = 2 cos(ω0 Ts )x(2) − x(1)
d(k) = x2 (k)

Como la perturbación d(k) no es conocida, siempre tendremos que usar un

estimador de estados para conocer la condición inicial que en este caso es:
· ¸
x(0)
X (0) =
ζ(0)

En un sistema aumentado cuyas matrices son:

· ¸ · ¸
A 0 B £ ¤
Aa = ; Ba = ; Ca = C Cζ .
0 Aζ 0