2.7: Ejercicios

Ejercicio 2.1 Ajuste de mínimos cuadrados de una elipse

Supongamos que un objeto en particular se modela como moviéndose en una órbita elíptica centrada en el origen. Su trayectoria nominal se describe en coordenadas rectangulares\((r, s)\) por la ecuación de restricción\(x_{1} r^{2}+ x_{2} s^{2}+ x_{3} rs=1\), donde\(x_{1}\)\(x_{2}\), y\(x_{3}\) son parámetros desconocidos que especifican la órbita. Tenemos disponibles las siguientes mediciones ruidosas de las coordenadas del objeto\((r, s)\) en diez puntos diferentes de su órbita:

\ [\ begin {array} {l}
(0.6728,0.0589) (0.3380,0.4093) (0.2510,0.3559) (-0.0684,0.5449)\\
(-0.4329,0.3657) (-0.6921,0.0252) (-0.3681, -0.2020) (0.0019, -0.3769)\\
(0.0825, -0.3508) (0.5294, -0.2918)
\ end {array}\ nonumber\]

Contestar

Se cree que las diez mediciones son igualmente confiables. Para tu comodidad, estos diez pares de\((r, s)\) valores medidos se han almacenado en vectores de columna nombrados\(r\) y a los\(s\) que puedes acceder a través del casillero 6.241 en Atena.* Después de agregar 6.241, y una vez en el directorio en el que estás ejecutando Matlab, puedes copiar los datos usando cp /MIT/6.241/Público/Fall95/HW1RS.MAT hw1rs.mat. Entonces, en Matlab, escriba load hw1rs para cargar los datos deseados; escriba quién para confirmar que los vectores\(r\) y de hecho\(s\) están disponibles.

Usando la ecuación de restricción asumida, podemos organizar la información dada en forma del sistema lineal de ecuaciones (aproximadas)\(A x \approx b\), donde\(A\) es una\(10 \times 3\) matriz conocida,\(b\) es un\(10 \times 1\) vector conocido, y\(x=\left(x_{1}, x_{2}, x_{3}\right)^{T}\). Este sistema de 10 ecuaciones en 3 incógnitas es inconsistente. Deseamos encontrar la solución\(x\) que minimice la norma euclidiana (o longitud) del error\(Ax - b\). Compare las soluciones obtenidas usando las siguientes cuatro invocaciones de Matlab, cada una de las cuales, en principio, da la solución de menor error cuadrado deseada:

a)\(x=A\backslash b \)
b) c\(x=\operatorname{pinv}(A) * b \)
) d\(x=\operatorname{inv}\left(A^{\prime} * A\right) * A^{\prime} * b\)
)\([q, r]=q r(A)\), seguida de la aplicación del enfoque descrito en el Ejercicio 3.1

Para obtener más información sobre estos comandos, prueba help slash, help qr, help pinv, help inv, etc. [Por cierto, el prime\(^{\prime}\),, en Matlab toma la transposición del complejo conjugado de una matriz; si quieres la transposición ordinaria de una matriz compleja\(C\), tienes que escribir \(C^{\prime}\)o\(transp(C)\).]

Debes incluir en tus soluciones una gráfica la elipse que corresponda a tu estimación de\(x\). Si crea el siguiente archivo de función en su directorio Matlab, con el nombre ellipse.m, puede obtener las coordenadas polares theta, rho de\(n\) puntos en la elipse especificada por el vector de parámetros\(x\). Para ello, ingrese [theta, rho] =elipse (x, n); en el prompt de Matlab. Luego puede trazar la elipse usando el comando polar (theta, rho).

\ [\ begin {array} {l}
\ text {función [theta, rho] =elipse (x, n)}\\
\%\\ texto {[theta, rho] = elipse (x, n)}\
\\%\
\%\ texto {El vector}\ x= [x (1), x (2), x (3)] ^ {\ prime}\,\ text {, define una elipse centrada en el origen}\\
\%\ text {vía la ecuación x (1) *}\ mathrm {r} ^ {\ wedge} 2 + x (2) *\ mathrm {s} ^ {\ wedge} 2+ x (3) *r*s=1\ text {.} \\%
\\ text {Esta rutina genera las coordenadas polares de puntos en el eclipse,}\\
\%\ text {para enviar a una orden de trazado. Lo hace resolviendo para el radial}\\
\%\ texto {distancia en n direcciones angulares igualmente espaciadas.} \\
\%\ text {Usa polar (theta, rho) para trazar realmente la elipse.} \\
\ end {array}\ nonumber\]
\ [\ begin {array} {l}
\ text {theta} =0:\ left (2^ {*}\ mathrm {pi}/\ mathrm {n}\ derecha):\ left (2^ {*}\ mathrm {pi}\ right);\
\\ mathrm {a} =\ mathrm {x} (1) {*}\ cos (\ texto {theta})\ cdot^ {\ cuña} 2+\ mathrm {x} (2) ^ {*}\ sin (\ text {theta})\ cdot^ {\ wedge} 2+\ mathrm {x} (3) ^ {*}\ izquierda (\ cos (\ text {theta})\ cdot^ {*}\ sin (\ texto {theta})\ derecha);\\
\ texto {rho} =\ nombreoperador {unos} (\ nombreoperador {tamaño} (\ mathrm {a}))\ cdot/\ mathrm {sqrt} (\ mathrm {a});
\ end {array}\ nonumber\]

Ejercicio 2.2 Aproximación por un polinomio

Vamos\(f(t)=0.5 e^{0.8 t}, t \in[0,2]\).

(a) Supongamos que 16 medidas exactas de\(f(t)\) están disponibles para usted, tomadas en los horarios\(t_{i}\) enumerados en la matriz T a continuación:

\ [\ izquierda. \ begin {array} {llllllll}
T= & {\ left [2\ cdot 10^ {-3},\ right.} & 0.136, & 0.268, & 0.402, & 0.536, & 0.668, & 0.802, & 0.936\\
& 1.068, & 1.202, & 1.336, & 1.468, & 1.602, & 1.736, & 1.868, & 2.000
\ end {array}\ derecho]\ nonumber\]

Utilice Matlab para generar estas mediciones:

\[y_{i}=f\left(t_{i}\right) \quad i=1, \ldots, 16 \quad t_{i} \in T\nonumber\]

Ahora determina los coeficientes de la aproximación polinomio de error de mínimos cuadrados de las mediciones, para

1. un polinomio de grado 15,\(p_{15}(t)\);
2. un polinomio de grado 2,\(p_{2}(t)\).

Comparar la calidad de las dos aproximaciones por trazado\(y(t_{i})\),\(p_{15}(t_{i})\) y\(p_{2}(t_{i})\) para todos\(t_{i}\) en T. Para ver qué tan bien estamos aproximando la función en todo el intervalo, también plot\(f(t)\),\(p_{15}(t)\) y\(p_{2}(t)\) en el intervalo [0, 2]. (Elija una cuadrícula muy fina para el intervalo, por ejemplo t= [0:1000] '/500.) Informe sus observaciones y comentarios.

(b) Ahora suponga que sus mediciones se ven afectadas por algún ruido. Generar las mediciones usando

\[y_{i}=f\left(t_{i}\right) + e(t_{i})\quad i=1, \ldots, 16 \quad t_{i} \in T\nonumber\]

donde el vector de valores de ruido se puede generar de la siguiente manera:

\ [\ begin {array} {l}
\ text {randn}\ left (^ {\ prime}\ text {seed} ^ {\ prime}, 0\ prime}, 0\ derecha);\\
e=\ nombreoperador {randn} (\ nombreoperador {siz} e (T));
\ end {array}\\ nonumber\]

Nuevamente determinar los coeficientes de la aproximación polinomio de error cuadrático mínimo de las mediciones para

1. un polinomio de grado 15,\(p_{15}(t)\);
2. un polinomio de grado 2,\(p_{2}(t)\).

Compare las dos aproximaciones como en la parte (a). Informe sus observaciones y comentarios. Explique cualquier resultado sorprendente.

(c) Hasta el momento hemos obtenido aproximaciones polinómicas de\(f(t), t \in [0, 2]\), aproximando las mediciones en\(t_{i} \in {T}\). Ahora nos interesa minimizar el error cuadrado de la aproximación polinómica en todo el intervalo [0, 2]:

\[\min \left\|f(t)-p_{n}(t)\right\|_{2}^{2}=\min \int_{0}^{2}\left|f(t)-p_{n}(t)\right|^{2} d t\nonumber\]

donde\({p}_{n}(t)\) hay algún polinomio de grado\(n\). Encuentra el polinomio\({p}_{2}(t)\) de grado 2 que soluciona el problema anterior. ¿Son muy diferentes entre sí lo óptimo\({p}_{2}(t)\) en este caso y lo óptimo\({p}_{2}(t)\) de las partes (a) y (b)? Elaborar.

Ejercicio 2.3 Combinar Estimaciones

Supongamos\(y_{1}=C_{1} x+e_{1}\) y\(y_{1}=C_{1} x+e_{1}\), donde x es un vector n, y\(C_{1}\),\(C_{2}\) tienen rango de columna completo. Dejar\(\widehat{x}_{1}\) denotar el valor de\(x\) que minimiza\(e_{1}^{T} S_{1} e_{1}\), y\(\widehat{x}_{2}\) denotar el valor que minimiza\(e_{2}^{T} S_{2} e_{2}\), donde\(S_{1}\) y\(S_{2}\) son matrices definitivas positivas. Demostrar que el valor\(\widehat{x}\) de\(x\) que minimiza se\(e_{1}^{T} S_{1} e_{1}+ e_{2}^{T} S_{2} e_{2}\) puede escribir enteramente en términos de\(\widehat{x}_{1}\),\(\widehat{x}_{2}\), y las\(n \times n\) matrices\(Q_{1}=C_{1}^{T} S_{1} C_{1}\) y\(Q_{2}=C_{2}^{T} S_{2} C_{2}\). ¿Cuál es el significado de este resultado?

Ejercicio 2.4 Estimaciones de ventana exponencialmente

Supongamos que observamos las medidas escalares

\[y_{i}=c_{i} x+e_{i}, \quad i=1,2, \ldots\nonumber\]

donde\(c_{i}\) y\(x\) son posiblemente vectores (vectores de fila y columna respectivamente).

(a) Mostrar (reduciendo esto a un problema que ya sabemos resolver - ¡no empieces de cero!) que el valor\(\widehat{x}_{k}\) de\(x\) eso minimiza el criterio

\[\sum_{i=1}^{k} f^{k-i} e_{i}^{2}, \quad \text { some fixed } f, \quad 0<f \leq 1\nonumber\]

está dado por

\[\hat{x}_{k}=\left(\sum_{i=1}^{k} f^{k-i} c_{i}^{T} c_{i}\right)^{-1}\left(\sum_{i=1}^{k} f^{k-i} c_{i}^{T} y_{i}\right)\nonumber\]

El llamado factor f de fundido u olvido nos permite ponderar preferencialmente las mediciones más recientes mediante picking\(0 < f < 1\), de manera que los datos antiguos se descuentan a una tasa exponencial. Entonces decimos que los datos han sido sometidos a desvanecimiento exponencial u olvido o ponderación o ventana o ahusamiento o.... Esto suele ser deseable, con el fin de mantener el filtro adaptativo a los cambios que puedan ocurrir en\(x\). De lo contrario el filtro se vuelve progresivamente menos atento a los nuevos datos y se queda dormido, con su ganancia acercándose a 0.

b) Demostrar ahora que

\[\hat{x}_{k}=\hat{x}_{k-1}+Q_{k}^{-1} c_{k}^{T}\left(y_{k}-c_{k} \hat{x}_{k-1}\right)\nonumber\]

donde

\[Q_{k}=f Q_{k-1}+c_{k}^{T} c_{k}, \quad Q_{0}=0\nonumber\]

El vector\(g_{k} = Q_{k}^{-1} c_{k}^{T}\) se denomina ganancia del estimador.

(c) Si\(x\) y\(c_{i}\) son escalares, y\(c_{i}\) es una constante\(c\), determinar\(g_{k}\) en función de\(k\). ¿Cuál es la ganancia en estado estacionario\(g_\infty\)? \(g_\infty\)Aumenta o disminuye a medida que\(f\) aumenta - y ¿por qué esperas esto?

Ejercicio 2.5

Supongamos que nuestro modelo para alguna forma de onda\(y(t)\)\(\alpha\) es\(y(t)=\alpha \sin (\omega t)\), donde es un escalar, y supongamos que tenemos mediciones\(y\left(t_{1}\right), \ldots, y\left(t_{p}\right)\). Debido a los errores de modelado y la presencia de ruido de medición, generalmente no encontraremos ninguna elección de parámetros del modelo que nos permita dar cuenta con precisión de todas las mediciones p.

(a) Si\(\omega\) se conoce, encuentre el valor de\(\alpha\) que minimice

\[\sum_{i=1}^{p}\left[y\left(t_{i}\right)-\alpha \sin \left(\omega t_{i}\right)\right]^{2}\nonumber\]

b) Determinar este valor de\(\alpha\) si\(\omega=2\) y si los valores medidos de\(y(t)\) son:

\ [\ begin {array} {llll}
y (1) =+2.31 & y (2) =-2.01 & y (3) =-1.33 & y (4) =+3.23\\
y (5) =-1.28 & y (6) =-1.66 & y (7) =+3.28 & y (8) =-0.88
\ end {array}\ nonumber\]

(Generé estos datos usando la ecuación\(y(t)=3 \sin (2 t)+ e(t)\) evaluada en los valores enteros\(t=1, \ldots, 8\), y\(t\) siendo\(e(t)\) para cada uno un número aleatorio distribuido uniformemente en el intervalo - 0.5 a +0.5.)

(c) Supongamos que\(\alpha\) y\(\omega\) son desconocidos, y que deseamos determinar los valores de estas dos variables que minimicen el criterio anterior. Supongamos que se le dan estimaciones iniciales\(\alpha_{0}\) y\(\omega_{0}\) para los valores de minimización de estas variables. Usando el algoritmo Gauss-Newton para este problema de mínimos cuadrados no lineales, es decir, aplicar LLSE al problema obtenido linealizando sobre las estimaciones iniciales, determinar explícitamente las estimaciones\(\alpha_{1}\) y\(\omega_{1}\) obtenidas después de una iteración de este algoritmo. Utilice la siguiente notación para ayudarle a escribir la solución en forma condensada:

\[a=\sum \sin ^{2}\left(\omega_{0} t_{i}\right), \quad b=\sum t_{i}^{2} \cos ^{2}\left(\omega_{0} t_{i}\right), \quad c=\sum t_{i}\left[\sin \left(w_{0} t_{i}\right)\right]\left[\cos \left(w_{0} t_{i}\right)\right]\nonumber\]

d) ¿Qué valores obtiene para\(\alpha_{1}\) y\(\omega_{1}\) con los datos dados en la letra b) anterior si las conjeturas iniciales son\(\alpha_{0}=3.2\) y\(\omega_{0}=1.8\)? Continuar la estimación iterativa unos pasos más. Repita el procedimiento cuando sean las conjeturas iniciales\(\alpha_{0}=3.5\) y\(\omega_{0}=2.5\), verificando que el algoritmo no converja.

(e) Dado que solo\(\omega\) ingresa al modelo de manera no lineal, podríamos pensar en un algoritmo descompuesto, en el que\(\alpha\) se estima usando mínimos cuadrados lineales y\(\omega\) se estima a través de mínimos cuadrados no lineales. Supongamos, por ejemplo, que nuestra estimación inicial de\(\omega\) es\(\omega_{0}=1.8\). Ahora obtenga una estimación\(\alpha_{1}\) de\(\alpha\) usar el método de mínimos cuadrados lineales que utilizó en (b). Después obtener una (¿mejorada?) estimación\(\omega_{1}\) de\(\omega\), usando una iteración de un algoritmo de Gauss-Newton (similar a lo que se necesita en (c), excepto que ahora solo está tratando de estimar\(\omega\)). A continuación, obtener la estimación\(\alpha_{2}\) a través de mínimos cuadrados lineales, y así sucesivamente. Compara tus resultados con lo que obtienes a través de este procedimiento descompuesto cuando tu estimación inicial sea\(\omega_{0}=2.5\) en lugar de 1.8.

Ejercicio 2.6 Comparando Estimadores Diferentes

Este problema le pide comparar el comportamiento de diferentes algoritmos de estimación de parámetros ajustando un modelo del tipo\(y(t)=a \sin (2 \pi t)+b \cos (4 \pi t)\) a datos ruidosos tomados a valores de\(t\) que están separados .02 en el intervalo (0,2].

Contestar

Primero sintetiza los datos sobre los que probarás los algoritmos. A pesar de que sus algoritmos de estimación asumirán eso\(a\) y\(b\) son constantes, estamos interesados en ver cómo rastrean los cambios de parámetros también. En consecuencia, vamos\(a = 2\),\(b = 2\) para los primeros 50 puntos, y\(a = 1\),\(b = 3\) para los siguientes 50 puntos. Para obtener (aproximadamente) variables aleatorias normalmente distribuidas, utilizamos la función randn para producir variables con media 0 y varianza 1.

Una manera elegante de generar los datos en Matlab, explotando la facilidad de Matlab con vectores, es definir los vectores t 1 = 0:02:0:02:1.0 y t 2 = 1:02:0:02:2.0, luego establecer

\[y 1=2 * \sin (2 * \mathrm{pi} * t 1)+2 * \cos (4 * \mathrm{pi} * t 1)+s * \operatorname{randn}(\operatorname{siz} e(t 1))\nonumber\]

y

\[y 2=\sin (2 * p i * t 2)+3 * \cos (4 * p i * t 2)+s * \operatorname{randn}(\operatorname{siz} e(t 2))\nonumber\]

donde s determina la desviación estándar del ruido. Elige\(s = 1\) para este problema. Por último, establecer\(y = [y1, y2]\). Sin bucles, sin contadores, sin alboroto!!

Ahora estime a y b a partir de y usando los siguientes algoritmos. Asumir estimaciones previas\(\widehat{a}_{0}= 3\) y\(\widehat{b}_{0}= 1\), ponderadas por igual con las mediciones (así todos los pesos se pueden tomar como 1 sin pérdida de generalidad). Trazar sus resultados para ayudar a la comparación.

i) Mínimos cuadrados recursivos

(ii) Mínimos cuadrados recursivos con memoria de desvanecimiento exponencial, como en el Problema 3. Uso\(f = .96\)

(iii) El algoritmo en (ii), pero con\(Q_{k}\) del Problema 3 sustituido por\(q_{k} = (1/n) \times trace(Q_{k})\), donde\(n\) está el número de parámetros, así que\(n = 2\) en este caso. (Recordemos que la traza de una matriz es la suma de sus elementos diagonales. Tenga en cuenta que\(q_{k}\) en sí mismo satisface una recursión, que debe anotar.)

(iv) Un algoritmo de la forma

\[\hat{x}_{k}=\hat{x}_{k-1}+\frac{.04}{c_{k} c_{k}^{T}} c_{k}^{T}\left(y_{k}-c_{k} \hat{x}_{k-1}\right)\nonumber\]

donde\(c_{k}=[\sin (2 \pi t), \cos (4 \pi t)]\) se evaluaron en el k ésimo instante de muestreo, así\(t = .02k\).

Ejercicio 2.7 Estimación recursiva de un vector de estado

Este curso pronto comenzará a considerar modelos estado-espacio de la forma

\[x_{l}=A x_{l-1}\ \tag{2.4}\]

donde\(x_{l}\) es un vector n que denota el estado en el momento\(l\) de nuestro modelo de algún sistema, y A es una\(n \times n\) matriz conocida. Por ejemplo, supongamos que el sistema de interés es una máquina giratoria, con posición angular\(d_{l}\) y velocidad angular\(\omega_{l}\) en el tiempo\(t = l T\), donde\(T\) hay algún intervalo de muestreo fijo. Si creyéramos que la máquina giraba a velocidad constante, seríamos conducidos al modelo

\ [\ left (\ begin {array} {l}
d_ {l}\
\ omega_ {l}
\ end {array}\ right) =\ left (\ begin {array} {ll}
1 & T\\
0 & 1
\ end {array}\ derecha)\ left (\ begin {array} {l}
d_ {l-1}\
\ omega_ {l-1}
\ end {array}\ derecha)\ nonumber\]

Supongamos que A es no singular a lo largo de este problema.

Para el ejemplo de máquina giratoria anterior, a menudo es de interés obtener estimaciones de menor error cuadrado de la posición y la velocidad (constante), utilizando mediciones ruidosas de la posición angular\(d_{j}\) en los instantes de muestreo. De manera más general, es de interés obtener una estimación de menor error cuadrado del vector de estado\(x_{i}\) en el modelo (2.4) a partir de mediciones ruidosas de componentes p-\(y_{j}\) que están relacionadas\(x_{j}\) con una ecuación lineal de la forma

\[y_{j}=C x_{j}+e_{j}, \quad j=1, \ldots, i\nonumber\]

donde C es una\(p \times n\) matriz. También asumiremos que se dispone\(\widehat{x}_{0}\) de\(x_{0}\) una estimación previa de:

\[\widehat{x}_{0}= x_{0}+ e_{0}\nonumber\]

Dejar\(\widehat{x}_{i|i}\) denotar el valor de\(x_{i}\) que minimiza

\[\sum_{j=0}^{i}\left\|e_{j}\right\|^{2}\nonumber\]

Esta es la estimación de\(x_{i}\) dada la estimación previa y mediciones hasta el tiempo\(i\), o la “estimación filtrada” de\(x_{i}\). Del mismo modo, vamos a\(\widehat{x}_{i|i-1}\) denotar el valor de\(x_{i}\) que minimiza

\[\sum_{j=0}^{i-1}\left\|e_{j}\right\|^{2}\nonumber\]

Esta es la estimación de menor error cuadrado de\(x_{i}\) dada la estimación previa y las mediciones hasta el tiempo\(i - 1\), y se denomina la “predicción de un paso” de\(x_{i}\).

a) Establecer el sistema lineal de ecuaciones cuya solución de error de mínimos cuadrados sería\(\widehat{x}_{i|i}\). De igual manera, configurar el sistema lineal de ecuaciones cuya solución de error de mínimos cuadrados sería\(\widehat{x}_{i|i-1}\).

b) Demostrar que\(\widehat{x}_{i|i-1}=A\widehat{x}_{i-1|i-1}\)

c) Determinar una recursión que se\(\widehat{x}_{i|i}\) exprese en términos de\(\widehat{x}_{i-1|i-1}\) y\(y_{i}\). Este es el prototipo de lo que se conoce como el filtro Kalman. Una versión más elaborada del filtro Kalman incluiría ruido aditivo que impulsa el modelo estado-espacio, y otros adornos, todo en un contexto estocástico (más que el determinista que se da aquí).

Ejercicio 2.8

Dejar\(\widehat{x}\) denotar el valor de\(x\) que minimiza\(\|y-A x\|^{2}\), donde\(A\) tiene rango de columna completa. Dejar\(\bar{x}\) denotar el valor de\(x\) que minimiza este mismo criterio, pero ahora sujeto a la restricción que\(z = Dx\), donde D tiene rango de fila completa. Demostrar que

\[\bar{x}=\hat{x}+\left(A^{T} A\right)^{-1} D^{T}\left(D\left(A^{T} A\right)^{-1} D^{T}\right)^{-1}(z-D \hat{x})\nonumber\]

(Sugerencia: Un enfoque para resolver esto es usar nuestra formulación recursiva de mínimos cuadrados, pero modificada para el caso limitante donde se sabe que uno de los conjuntos de medición, es decir,\(z = Dx\) en este caso, no tiene error. Es posible que tenga que usar algunas de las identidades matriciales del capítulo anterior).