Un sistema de ecuaciones matriciales (2 Riccati, 1 Lyapunov)

Question

Preparado:

Dejemos que $\gamma \in(0,1)$ , ${\bf F},{\bf Q} \in \mathbb R^{n\times n}$ , ${\bf H}\in \mathbb R^{n\times r}$ y ${\bf R}\in \mathbb R^{r\times r}$ y supongamos que ${\bf P}$ , ${\bf W}$ , ${\bf X}\in \mathbb R^{n\times n}$ y ${\bf K}$ , ${\bf L}\in \mathbb R^{n\times r}$ satisfacer

\begin{align} {\bf P} &={\bf F}({\bf I}_{n}-{\bf K} {\bf H}^\top){\bf P}{\bf F}^\top+{\bf Q}, \;\;\;\;\;\;\;\;\text{where}\;\;\;\;{\bf K}\equiv {\bf P} {\bf H}\left({\bf H}^\top{\bf P} {\bf H}+ \frac{1}{\gamma}{\bf R} \right)^{-1} \tag1\\[4ex] {\bf W} &={\bf F}({\bf I}_{n}-{\bf L} {\bf H}^\top){\bf W}{\bf F}^\top+{\bf Q}, \;\;\;\;\;\;\;\;\text{where}\;\;\;\;{\bf L}\equiv {\bf W} {\bf H}({\bf H}^\top {\bf W} {\bf H}+ {\bf R})^{-1} \tag2\\[4ex] {\bf X} &={\bf K}{\bf H}^\top {\bf W}+({\bf I}_n-{\bf K}{\bf H}^\top){\bf F}{\bf X}({\bf I}_n-{\bf H} {\bf L}^\top){\bf F}^\top \tag3\\ {\color{white}X} \end{align}

Además, supongamos que ${\bf P}$ , ${\bf W}$ , ${\bf R}$ y ${\bf Q}$ son simétricos, y ${\bf P}$ , ${\bf W}$ y ${\bf X}$ son invertibles.

---

Quiero probarlo:

$${\bf K}=(\gamma{\bf I}_{n}+(1-\gamma) {\bf X} {\bf W}^{-1}){\bf L} \tag4$$

---

Algunas ideas y comentarios:

• En el caso escalar (con $n=r=1$ ) lo que funciona es restar $(2)$ de $(1)$ que elimina ${\bf Q}$ y luego utilizarlo para resolver ${\bf F}^2$ es decir, (suprimiendo el $^\top$ notación) $$ {\bf F}^2 = \frac{(\gamma{\bf H}^2{\bf P}+{\bf R})({\bf H}^2{\bf W}+{\bf R})({\bf P}-{\bf W})}{{\bf R}^2({\bf P}-{\bf W})+(1-\gamma){\bf H}^2{\bf R}{\bf P}{\bf W})},$$ sustituir esto por $(3)$ y resolver para ${\bf X}$ que da como resultado $${\bf X} = \frac{(\gamma {\bf R}({\bf P}-{\bf W})+\gamma(1-\gamma){\bf H}^2{\bf P}{\bf W}}{(1-\gamma)(\gamma {\bf H}^2{\bf P}+{\bf R}))}. $$ Finalmente, resolviendo para ${\bf X}$ utilizando $(4)$ produce lo mismo.

• Obsérvese que las ecuaciones $(1)$ y $(2)$ son ecuaciones de Riccati en ${\bf P}$ y ${\bf W}$ respectivamente, lo que significa que existen métodos conocidos (descritos aquí ) para resolver ${\bf P}$ y ${\bf W}$ . Aunque no pude hacer uso de esas soluciones.

• Ecuación $(3)$ es similar a una ecuación de Lyapunov en ${\bf X}$ lo que significa que un método de vectorización (descrito aquí ) permite obtener una ecuación para ${\bf X}$ .

• Creo que la prueba vendrá de un procedimiento similar al que funciona para el caso escalar. Es decir, restando $(2)$ de $(1)$ utilizando esta ecuación para eliminar ${\bf F}$ de $(3)$ y luego mostrar que la ecuación resultante implica $(4)$ .

• Si puedes probarlo bajo supuestos adicionales, eso también podría ser útil.

---

A continuación se presenta un sencillo código de Matlab que permite comprobar el resultado:

    gamma = 0.5;

    F = [2, 2, 3; 4, 5, 6; 7, 8, 9];  % n x n
    H = [1, 2; 3, 1; 2, 1];           % n x r
    R = [3 , 1; 1, 3];                % r x r, symmetric
    Q = [2, 0, 0; 0, 5, 0; 0, 0, 9];  % n x n, symmetric

    n = length(F); 
    r = length(R); 

    I = eye(n);

    P = I;
    W = I;
    X = I;
    for i=1:1000
        K = (P*H)/(H'*P*H+(R/gamma));
        P = F*(P-K*H'*P)*F'+Q;
        L = (W*H)/(H'*W*H+R);    
        W = F*(W-L*H'*W)*F'+Q;
        X = K*H'*W+(I-K*H')*F*X*(I-H*L')*F';
    end

    disp(['K - (gamma*I+(1-gamma)*X*W^(-1))*L = ',...
          num2str(sum(sum(abs(K-(gamma*I+(1-gamma)*X*(W^(-1)))*L))))])

---

Esta pregunta es una versión reducida de otra pregunta: Equivalencia de previsión entre dos filtros de Kalman de estado estacionario (con ${\bf L} \equiv {\bf L}_1$ , ${\bf W} \equiv {\bf W}_{11}$ y ${\bf X} \equiv {\bf W}_{12}$ ).

Answer 1

Lo pongo como otra respuesta porque no quiero borrar nada de lo anterior todavía y la otra opción es hacerlo demasiado largo.

4) La ecuación para $X$ es Ecuación de Sylvester pero tienes razón esto está estrechamente relacionado con Lyapunov. La ecuación de Lyapunov es un caso especial de la ecuación de Sylvester.

5) Eliminar la etiqueta de álgebra lineal, el problema es no lineal.

6) Añade etiquetas de optimización y control óptimo.

7) Ya que hemos averiguado que el problema de tu código es la conocida inestabilidad numérica de la inversa de la matriz (debería haberlo notado al instante). He tratado de escribir un código que evite el problema, es decir, que no invierta las matrices explícitamente. He terminado con el siguiente código que también da más pistas sobre dónde buscar condiciones adicionales sobre las matrices. Esperemos que estas condiciones ayuden a demostrar la conjetura.

Brevemente sobre el algoritmo,

• utilizar el algoritmo de matlab dare para obtener las matrices de ricatti $P$ y $W$
• Para obtener los valores de $K$ y $L$ utilizar el algoritmo de búsqueda del cero ( fsolve ) utilizando su fórmula en ( $1$ ) y ( $2$ ) respectivamente.
• Utilizar el algoritmo de Matlab dlyap para obtener $X$
• Para verificar el resultado resolver perturbar la ecuación requerida para que no haya una matriz inversa, es decir $${\bf K}=(\gamma{\bf I}_{n}+(1-\gamma) {\bf X} {\bf W}^{-1}){\bf L}$$ $${\bf K}=(\gamma{\bf I}_{n}+(1-\gamma) {\bf X} {\bf W}^{-1}) {\bf W} {\bf H}({\bf H}^\top {\bf W} {\bf H}+ {\bf R})^{-1} $$ $${\bf K}({\bf H}^\top {\bf W} {\bf H}+ {\bf R})=(\gamma{\bf I}_{n}+(1-\gamma) {\bf X} {\bf W}^{-1}) {\bf W} {\bf H} $$ $${\bf K}({\bf H}^\top {\bf W} {\bf H}+ {\bf R})=(\gamma{\bf W}+(1-\gamma) {\bf X} ) {\bf H} $$

---

posibles pistas para la conjetura El algoritmo es más estable con respecto a la inversa de la matriz. El hallazgo de $K,L$ es un poco un problema pensado porque la conjetura inicial no está guiada. En general, puede ser esencial la dificultad para llegar a los valores correctos. Sin embargo, cuando ejecute el algoritmo $1000$ veces que me $>600$ falla de dare y sólo unos pocos fallos de fsolve . Creo que dare puede arrojar algo de luz sobre el problema original.

Así, he descubierto que dare es muy caprichoso algoritmo. Más concretamente, no tiene una solución para cualquier entrada, sino para una familia reducida de matrices. Las condiciones exactas como se describe en matlab no son claras para mí, por lo que incluso pidió otra pregunta con la esperanza de que alguien dé una respuesta. He encontrado varias preguntas sin respuesta sobre las ecuaciones algebraicas de Ricatti, por lo que sugeriría acudir a la literatura para descubrir en qué condiciones $F,Q,H,R$ la condición de que la solución de Ricatti existe. Lo mismo sobre Silverster\Lyapunov. También investigaría la relación entre las ecuaciones algebraicas y las diferenciales, tal vez la reducción a DE facilite la demostración, así que vale la pena intentarlo.

---

El código

n=9;
r=5;

F=round(10*rand(n),0);
H=round(10*rand(n,r),0);
R= round(10*rand(r),0);
R=R+R';
Q=round(10*rand(n),0);
Q=Q+Q';

P = dare(F',H,Q,R/gamma);
W = dare(F',H,Q,R);

options=optimset('Algorithm',{'levenbergmarquardt',0.0005});%,'Display','iter');
[K,fvalK,exitflagK] = fsolve(@(x) F*(P-x*H'*P)*F'+Q-P,zeros(n,r),options);
[L,fvalL,exitflagL] = fsolve(@(x) F*(W-x*H'*W)*F'+Q-W,zeros(n,r),options);

if exitflagK<=0 || exitflagL<=0 
        error('try again') 
end

X = dlyap((I-K*H')*F,(I-H*L')*F',K*H'*W);

tmp = K*(H'*W*H+R)-(gamma*W+(1-gamma)*X)*H;   
err=norm(tmp(:),inf);
disp(['K(H''W*H+R) - (gamma*W+(1-gamma)*X)*H) = ', num2str(err)])

Answer 2

1) Espero que me equivoque, pero he hecho un pequeño cambio en su programa (aleatorizado las matrices $F,H,R,Q$ ) y obtuve un resultado muy interesante: tu conjetura es falsa. Supongo que hay restricciones adicionales que te faltan, pero no tuve éxito en averiguarlas. Si ninguna de las dos cosas es cierta - por favor, mira los cambios que hice y hazme saber lo que estaba mal. Tenga en cuenta que con pequeñas $n$ y $r$ los valores tienden a ser pequeños para la mayoría de las ejecuciones, pero para las matrices grandes cada ejecución resulta en un número grande.

2) Usted dice que $X$ es similar a la matriz de Lyapunov, pero no logré reconocerla. En particular, hay una Matriz hermitiana ¿cuál de ellos se supone que es hermitiano? ¿Podría mostrar esta similitud explícitamente?

3) Intenté utilizar las identidades que puedes encontrar en los siguientes enlaces, pero no llegué muy lejos, pero tampoco le dediqué el tiempo suficiente.

Inversa de la suma de matrices

https://ecommons.cornell.edu/bitstream/handle/1813/32750/BU-647-M.version2.pdf?sequence=1

gamma = 0.5;
n=15;
r=7;

F=round(100*rand(n),0);  % n x n
H=round(100*rand(n,r),0);           % n x r
R= round(100*rand(r),0); %r x r,                 
R=R+R'; % make R symmetric

Q=round(100*rand(n),0);  % n x n, 
Q=Q+Q';% make Qsymmetric

I = eye(n);

P = I;
W = I;
X = I;
for i=1:1000
        %K = P*H*(H'*P*H+(R/(1-alpha)))^(-1);
        K = P*H*(H'*P*H+(R/(gamma)))^(-1);
        P = F*(P-K*H'*P)*F'+Q;
        L = W*H*(H'*W*H+R)^(-1);    
        W = F*(W-L*H'*W)*F'+Q;
        X = K*H'*W+(I-K*H')*F*X*(I-H*L')*F';
end

    disp(['K - (gamma*I+(1-gamma)*X*W^(-1))*L = ',...
          num2str(sum(sum(abs(K-(gamma*I+(1-gamma)*X*(W^(-1)))*L))))])
  F
  H
  R
  Q

El resultado de una carrera con $n=5$ , $r=2$

K - (gamma*I+(1-gamma)*X*W^(-1))*L = 2612.9753

F =

    80    64    11    68    56
    41    59    12    83    83
    97    67    35    12    49
    98    47    80    55    46
    64    12     2    78    56

H =

    64    52
    54    21
    98    54
    94    49
    46    28

R =

    24   124
   124    38

Q =

   164   106    60    77   105
   106    44   102    64   119
    60   102   138   108    38
    77    64   108   150    65
   105   119    38    65   104