Diferencia entre revisiones de «Aproximación de problemas de control.Derivadas topológicas para la resolución problemas inversos: Resolución actividad propuesta»

Revisión actual del 10:22 25 mar 2015

Plantilla:Artículo

El texto que aparece a continuación esta basado en la presentación Derivadas topológicas para la resolución de problemas inversos realizada por M. L. Rapún en la Escuela Técnica Superior de Ingenieros de Caminos el 11 de Febrero de 2015.

1 Introducción

1.1 General

Un problema inverso es aquel en el que a partir de ciertas mediciones o valores observados se debe deducir alguno o varios de los parámetros del sistema.

En este ejercicio, vamos a hacer énfasis en aquel problema inverso en el cual el parámetro que se pretende estimar son los defectos o inhomogeneidades del medio. Concretamente, se trataría de un problema de dispersión, en el cual se hace incidir una onda [math]u_{inc}[/math] en el seno de un medio [math]R[/math] contra unos defectos [math]\Omega[/math]. Dado que las propiedades del medio difieren en el defecto, el dominio [math]\Omega[/math] puede ser estimado por un procedimiento más o menos complejo.

El objetivo en este problema inverso concreto, es, a partir de las mediciones [math]u_{meas}[/math] tomados en una serie de receptores [math]\Gamma[/math], encontrar los objetos [math]\Omega[/math] tales que la diferencia [math]|u-u_{meas}|[/math] sea mínima.

Se debe subrayar que a pesar de su enorme interés, se trata de un problema mal planteado que puede no tener solución, y si la tiene puede no depender continuamente de los datos.

Este problema tiene numerosas aplicaciones, entre las que destacan aquellas en campos de tanto interés como puedan ser la medicina, la geología o el análisis de estructuras.

1.2 Problema modelo

Supongamos que estamos tratando con ondas acústicas. Tenemos que

[math]\rho U_{tt} = k \Delta U,[/math]

donde [math] k [/math] es la constante elásticas en problemas de dispersión acústica. La onda incidente es armónica en el tiempo

[math]U_{inc}(\vec{x}, t) = Re\left[e^{−i \omega t} u_{inc}(\vec{x})\right],[/math]

y [math]u_{inc}[/math] es una onda plana que se propaga en la dirección [math]\vec{d}[/math], tal que

[math]u_{inc}(\vec{x})= e^{i k \vec{x} \cdot \vec{d}}.[/math]

La solución del problema directo es también armónica en el tiempo por lo que

[math]U(\vec{x}, t) = Re\left[e^{−i \omega t} u(\vec{x})\right].[/math]

La onda incidente genera una la parte dispersada en [math]R \backslash \Omega[/math] y en una parte trasmitida en [math]\Omega[/math]

[math]u=u_{inc}+u_{sc} \text{ en } R\backslash\Omega[/math]

y

[math]u=u_{trans} \text{ en } \Omega.[/math]

Se tiene que $u$ es la solución del problema,

[math] \Delta u + k_e u = 0 \text{ en } R \backslash \Omega, [/math]

[math] \Delta u + k_i u = 0 \text{ en } \Omega, [/math]

con las condiciones de contorno

[math] u^+ = u^-\text{ en } \partial\Omega,[/math]

[math] \partial_n u^+ = \partial_n u^- \text{ en } \partial\Omega,[/math]

[math] \lim_{r \to \infty} r^{(n-1)/2} \left(\partial_r (u -u_{inc}) - i k_e (u -u_{inc})\right) = 0,[/math]

donde [math]k_e[/math] y [math]k_i[/math] corresponden al medio y a los defectos respectivamente, [math]r[/math] es el radio.

2 Derivadas topológicas

2.1 General

La variable que pretendemos encontrar, [math]\Omega[/math], puede ser hallada minimizando un funcional del error

[math] J=\dfrac{1}{2} \int_\Gamma |u-u_{meas}|^2.[/math]

La minimización puede ser realizada utilizando, entre otros, el método de las Derivadas topológicas. La Derivada topológica se define como,

[math]D_T(\vec{x}, R)=\dfrac{J(R\backslash B_\epsilon(\vec{x}))-J(R)}{f(\epsilon)},[/math]

donde [math]B_\epsilon(\vec{x})[/math] es la bola de centro [math]\vec{x}[/math] y radio [math]\epsilon[/math] y [math]f(\epsilon)[/math] es una función positiva monótona decreciente cuyo límite existe, es no nulo y cuyo valor tiende a cero cuando [math] \epsilon [/math] tiende a cero.

Ciertamente, la Derivada topológica es una función lineal de [math]\vec{x}[/math] fácil de interpretar como la variación del funcional [math]J[/math] cuando se retira un volumen infinitesimal centrado [math]\vec{x}[/math] y de radio [math]\epsilon[/math]. En forma de desarrollo,

[math] J(R\backslash B_\epsilon(\vec{x}))=J(R)+D_T(\vec{x}, R)f(\epsilon)+\ldots [/math]

Una formulación que permite entender mas claramente que en el caso en el que [math]D_T(\vec{x}, R)\lt0 [/math] la probabilidad de que exista un defecto localizado en [math]\vec{x}[/math] es elevada.

2.2 Calculo de la derivada topologica

En la práctica no se calcula la [math]D_T[/math] utilizando su definición sino que se utilizan formulaciones alternativas utilizando, por ejemplo, una definición equivalente como límite de derivadas de forma y realizando posteriormente desarrollos asintóticos. En la mayoría de los casos, se suele poder llegar a expresiones sencillas de la DT que involucran la resolución de un problema directo y de un problema adjunto. Esto se puede concretizar de la siguiente forma. Para

[math] J=\dfrac{1}{2} \int_\Gamma |u-u_{meas}|^2,[/math]

se tiene que

[math] D_T(\vec{x}, R)=Re\left[(k_i^2-k_e^2) u(\vec{x})w(\vec{x})\right] [/math]

donde [math]u[/math] es solución de un problema directo y [math]w[/math] es solución de un problema adjunto. Es decir, tomando [math]\Omega=\emptyset[/math] se tiene que [math]u[/math] es solución de

[math] \Delta u + k_e u = 0 [/math]

y por tanto es

[math]u_{inc}(\vec{x})= e^{i k \vec{x} \cdot \vec{d}}.[/math]

Por otra parte, [math]w[/math] es solución de

[math] \Delta w+ k_e^2 w = (\overline{u-u_{meas}})\partial_{\Gamma}.[/math]

Por ello, [math]w[/math] es del tipo

[math] w = \int_{\Gamma} G(\vec{x}-\vec{y}) (\overline{u-u_{meas}})(\vec{y}) dl_y [/math] donde [math]G(\vec{x}-\vec{y})[/math] es el Kernel de una convolución.

3 Aplicación de las derivadas Topológicas

3.1 Descripción del problema propuesto

Como parte del curso de Ecuaciones en derivas parciales, Problema inverso y control, se pide aplicar la técnica anteriormente descrita a un problema.

En una primera fase, el problema propuesto consiste en encontrar la derivada topológica en un dominio cuadrado [math][-1,1] \times [-1,1][/math] en el que existen unos defectos desconocidos. Para ello se han hecho incidir diez ondas con diversas direcciones y se cuenta con una doce puntos de medida, localizados en torno al dominio cuadrado, ver diagrama adjunto.

Para la resolución del problema se nos facilita un fichero en el que están contenidos:

• Las coordenadas de los puntos de medida
• Las medidas en los transductores
• Los ángulos de las ondas incidentes
• El numero de onda de la ecuación de Helmholtz

Con objeto de estudiar la sensibilidad de la solución a diferentes posiciones de los transductores, se pide además del caso ya descrito los siguientes resultados adicionales:

• Obtener la Derivada topológica utilizando los datos recogidos por los seis primeros receptores
• Obtener la Derivada topológica utilizando los datos recogidos en los seis últimos receptores
• Obtener la Derivada topológica utilizando los datos recogidos por los tres primeros receptores
• Obtener la Derivada topológica utilizando los datos recogidos por uno de cada dos receptores
• Obtener la Derivada topológica utilizando los datos recogidos en uno de cada tres receptores

Debido a que en el diagrama en el que se propone el problema existe el caso en el que se tienen en cuenta tan solo tres receptores (uno de cada cuatro) se añade este caso adicionalmente.

3.2 Metodo de solucion simplificada

Para encontrar la Derivada topológica vamos a calcular la funciónes [math]u[/math] y [math]w[/math] soluciones del problema directo e inverso.

Las ondas planas incidentes pueden ser caracterizadas por las formulas [math] u_{inc \, k} (\vec{x})=e^{i \lambda \vec{d_k} \cdot \vec{x}} [/math]

donde por simplicidad tomamos los

[math] d_{k}=(\cos{(\theta_k)}, \sin{(\theta_k)}). [/math]

donde [math]k[/math] representa cada una de las direcciones de la onda incidente.

En este ejercicio simplificado se ha tomado un valor de [math]\lambda=4[/math].

El funcional [math]J[/math] que debemos minimizar toma la expresión,

[math] J=\dfrac{1}{2} \sum_{k=1}^{N_\text{direcciones}}\sum_{j=1}^{N_\text{receptores}} |u_k(x_j)-u_{meas}(x_j)|^2.[/math]

La función [math]u_k[/math] es solución de

[math] \Delta u_k + \lambda u_k = 0 \text{ en } R^2 \backslash \Omega [/math]

con

[math] u_k = 0 \text{ en } \partial\Omega[/math]

y

[math] \lim_{r \to \infty} r^{(n-1)/2} \left(\partial_r (u_k -u_{inc \, k}) - i k_e (u_k -u_{inc \, k})\right) = 0.[/math]

La Derivada topológica resulta entonces ser

[math] D_T (\vec{x}) = \sum_{k=1}^{N_\text{direcciones}} Re\left[u_{inc \, k}(\vec{x}) w_k(\vec{x})\right]. [/math]

La solución del problema inverso se puede tomar de la forma

[math] w_k (\vec{x}) = \sum_{j=1}^{N_\text{receptores}} \dfrac{i}{4} H_0^{(1)}(\lambda |\vec{x}-\vec{x}_j|)(\overline{u_{meas \, k}(\vec{x}_j)-u_{inc \, k}(\vec{x}_j)}),[/math]

donde [math]H_0^{(1)}[/math] es la función de Hankel de primera especie y orden cero.

3.3 Código utilizado

La Derivada topológica fue calculada con el Script escrito en el lenguaje de programación Python que se muestra a continuación. Con el objeto de clarificar el codigo se han insertado comentarios que clarifican la función de cada bloque de código.

#!/usr/bin/python

# Cargar los paquetes necesarios

import scipy.io
import numpy as np
import scipy.special as sci_spc
import math

#Leer el archivo con los datos del probelma

mat = scipy.io.loadmat('DatosProblema.mat')

#Copiar cada variable por separado para mayor claridad

Xreceptores=np.array(mat['Xreceptores'])
Yreceptores=np.array(mat['Yreceptores'])
Angulos=np.array(mat['Angulos'])
OndaMedida=np.matrix(mat['OndaMedida'])
ll=np.array(mat['lambda'])
lam=ll[0,0]

#Generar direcciones. tomo como numero complejo aunque se trata de R^2

d = np.zeros(shape=Angulos.shape, dtype=np.complex)
for i in np.arange(Angulos.shape[1]):
    d[0,i]=math.cos(Angulos[0,i])+math.sin(Angulos[0,i])*1j

#Generar matriz de puntos para calcular

x =  np.zeros([201,201], dtype=np.complex)
for i in np.arange(-100,100,1):
    for j in np.arange(-100,100,1):
        x[i,j]=np.float(i)/100+np.float(j)/100*1j

# Generar matriz de puntos para la Derivada topologica

DT =  np.zeros([201,201], dtype=np.float)

# Generar vector con las coordenadas de los trasnductores

xj =  np.zeros([Xreceptores.shape[1]], dtype=np.complex)
for i in np.arange(Xreceptores.shape[1]):
    xj[i]=Xreceptores[0, i]+1j*Yreceptores[0, i]

# Valor de u en los trasnductores para cadqa direccion y trasnductor
uj = np.zeros(shape=OndaMedida.shape, dtype=np.complex)
for i in np.arange(uj.shape[0]):
    for k in np.arange(uj.shape[1]):
        uj[i,k]=np.exp(1j*lam*(d[0,k].real*xj[i].real+d[0,k].imag*xj[i].imag))


#print x[-90,-90]

# Generar matrices para u y w

u =  np.zeros([201,201,Angulos.shape[1]], dtype=np.complex)
w =  np.zeros([201,201,Angulos.shape[1]], dtype=np.complex)
ww =  np.zeros([201,201,Angulos.shape[1],Xreceptores.shape[1]], dtype=np.complex)

#Calcular u y w en las coordenadas para cada direccion y receptor

for a in np.arange(-100,100,1):
    for b in np.arange(-100,100,1):
        for k in np.arange(Angulos.shape[1]):
            u[a,b,k]=np.exp(1j*lam*(d[0,k].real*x[a,b].real+d[0,k].imag*x[a,b].imag))
            for j in np.arange(Xreceptores.shape[1]):
                ww[a,b,k,j]=1j/4*sci_spc.hankel1(0, lam*np.absolute(x[a,b]-xj[j]))*np.conjugate(OndaMedida[j,k]-uj[j,k])

# Sumar en receptores

for a in np.arange(-100,100,1):
    for b in np.arange(-100,100,1):
        for k in np.arange(Angulos.shape[1]):
             for j in np.arange(Xreceptores.shape[1]):
                 w[a,b,k]+= ww[a,b,k,j]


# Sumar en cada direccion

for a in np.arange(-100,100,1):
    for b in np.arange(-100,100,1):
        for k in np.arange(Angulos.shape[1]):
            DT[a,b]+=(u[a,b,k]*w[a,b,k]).real

# Output de los valores

for a in np.arange(-100,100,1):
    for b in np.arange(-100,100,1):
        #print (%f, %f, %f) % (x[a,b].real(), x[a,b].imag(), DT[a,b])        
        print '%f, %f, %f' % (x[a,b].real, x[a,b].imag, DT[a,b])

Los valores obtenidos fueron dibujados en forma de Mapa térmico utilizando un segundo Script, que se muestra a continuación, también realizado en lenguaje Python

#!/usr/bin/python

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.mlab import griddata

xi = np.linspace(-1., 1., 200)
yi = np.linspace(-1., 1., 200)

do=np.loadtxt('datos_escalar.dat', delimiter=',')

x=do[:,0]
y=do[:,1]
z=do[:,2]

zi = griddata(x, y, z, xi, yi, interp='linear')
CS = plt.contour(xi, yi, zi, 15, linewidths=0.5, colors='k')
CS = plt.contourf(xi, yi, zi, 15, cmap=plt.cm.rainbow,
                  vmax=abs(zi).max(), vmin=-abs(zi).max())
cbar = plt.colorbar() 
#cbar.set_label('DT', rotation=270)
cbar.set_label('DT, [-]')

plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xlabel('x, [-]')
plt.ylabel('y, [-]')

plt.show()

4 Soluciones Obtenidas

4.1 Derivada topológica obtenida teniendo en cuenta la totalidad de receptores

Las soluciones obtenidas teniendo en cuenta todos los receptores, localizados según el diagrama, se pueden observar en el diagrama térmico siguiente. Dos zonas diferenciadas aparecen claramente en la figura en las que [math]D_T[/math] toma valores negativos. Estas áreas están caracterizadas por adoptar tonos azules.

El diagrama obtenido sirve de caso base para comparar los resultados que se obtiene cuando se utilicen un numero inferior de receptores o una disposición distinta (objeto de los siguientes subapartados).

4.2 Derivada topológica obtenida teniendo en cuenta 6 receptores agrupados

Seguidamente se realizó una repetición del calculo pero teniendo en cuenta tan solo la mitad de los receptores. El calculo se realizo tomando inicialmente la primera mitad de los receptores y luego la segunda mitad. Para lograr una mayor claridad en la exposición se incluyen sendos diagramas en los que se puede apreciar la distribución de los receptores en cada uno de los casos.

Los resultados obtenidos parecen, a primera vista, casi simétricos con respecto a los receptores eliminados.

El resultado para el defecto mas alejado de los transductores se degrada notablemente en ambos caso. Esto puede ser debido probablemente a la existencia de un cierto apantallamiento por parte del defecto mas cercano a los transductores.

Los niveles que se obtienen hacen pensar, teniendo en cuenta el resultado de la sección anterior, que el defecto (o inclusión) inferior derecho tiene un tamaño superior que la inclusión superior izquierda.

4.3 Derivada topológica obtenida teniendo en cuenta 3 receptores agrupados

El procedimiento fue repetido teniendo en cuenta tan solo tres transductores, localizados según aparece en la figura.

El resultado muestra que la solución obtenida se ha seguido degradando. Los defectos aparecen con una forma significativamente alargada, aunque, probablemente debido a la localización de los obstáculos, su localización toma, como en el caso base, una forma simétrica.

4.4 Derivada topológica obtenida teniendo en cuenta receptores no agrupados

Por ultimo se tratan los resultados obtenidos en los casos en los que se han eliminado uno de cada dos, uno de cada tres y uno de cada cuatro receptores. Para mayor claridad se incluyen los diagramas correspondientes a cada una de las disposiciones de los transductores. Es visible como la imagen obtenida se va degradando poco a poco, especialmente en el últimos caso, en el que los obstáculos dejan de ser redondeados para achantarse perceptiblemente.

Dominio de calculo (granate) puntos de medida (estrellas) y ondas incidentes (flechas) para el problema inverso que se pide resolver