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 del 13:36 17 mar 2015

Plantilla:Artículo

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 de un 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.

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

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

donde [math] k [/math] es la constantes 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]

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

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

donde [math]k_e[/math] y [math]k_i[/math] corresponden al medio y a los defectos respectivamente.

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 el método de las Derivadas topológicas.

La Derivada topológicas 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 y es no nulo y que además tiende a cero cuando [math] \epsilon [/math] también 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]. Escrito en forma de desarrollo,

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

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]

3 Aplicación de las derivadas Topológicas

3.1 Descripción del problema propuesto

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

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, ver diagrama adjunto. 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.

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 con la configuración ya descrita las siguientes soluciones adicionales:

• Utilizando los datos recogidos por los 6 primeros receptores
• Utilizando los datos recogidos en los 6 últimos receptores
• Utilizando los datos recogidos por los 3 primeros receptores
• Utilizando los datos recogidos por uno de cada dos receptores
• 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ón [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} 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. Además, 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]

y con

[math] u_k = 0 \text{ en } \partial\Omega.[/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]

y tomando la solución del problema inverso como

[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]

en esta ultima ecuación [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 un Script en el lenguaje de programación Python.

#!/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 termico con una segundo Script también realizado en Python

3.4 Soluciones

Res1

Res1

Res1

Res1

Res1

Res1

Res1