Diferencia entre revisiones de «Explotación minera (grupo5-A)»

Revisión del 00:43 6 mar 2015

INTRODUCCIÓN

En una determinada región se ha encontrado el yacimiento de un mineral altamente demandado por la población, por ello vamos a tratar de modelizar el comportamiento de la explotación y con los resultados obtenidos analizaremos la veracidad del método empleado. Un trabajo de modelización consiste en, a partir de unos datos obtenidos de la realidad, obtener una caracterización matemática de dichos datos para poder aproximar su posible comportamiento siempre teniendo en cuenta las hipótesis de partida y teniendo en cuenta que hay que ser críticos, como veremos con los resultados obtenidos.

DATOS E INFORMACIÓN REAL

Mediante prospecciones geológicas se estima la cantidad total de mineral extraíble en 10875 Toneladas--------Variable K.

Un estudio de mercado dictamina que se espera un crecimiento rápido de la producción hasta los primeros 25 años, momento a partir del cual, ésta empieza a descender. Por las características técnicas de la explotación se estima una producción máxima de 240 Toneladas/año.

Además la mínima producción para la que la explotación es rentable son 25T/año, a partir de este momento se considera que habrá terminado la vida útil de la explotación.

ELECCIÓN DEL MÉTODO DE MODELIZACIÓN DEL PROCESO

Debido a estas características se decide utilizar el modelo de Gompertz que corresponde a la ecuación diferencial siguiente:

                                                           dQ/dt=rQLn(k/Q)

Donde:

Q(t): Cantidad de mineral extraída hasta el tiempo t [Toneladas]

k:Cantidad de mineral extraible [Toneladas]

r:Tasa intrínseca de crecimiento (Factor a determinar)

Para que nuestra ecuación diferencial quede completamente definida calculamos el valor de la constante r,Como tenemos datos reales sobre la producción y no sobre la cantidad de mineral extraído definimos la función producción como la dQ/dt es decir, la cantidad de mineral que se extrae por unidad de tiempo, lo que se corresponde perfectamente con la definición de producción, así:

                                                             P(Q)=rQLn(k/Q)

De esta función es importante decir que no depende de la variable tiempo implicitamente sino atraves de la función Q(t), es decir, la función producción es una composición de la función Q(t), por lo que, como es lógico pensar, depende del tiempo aunque no de manera explícita.

CÁLCULO DE r PARA DEFINIR COMPLETAMENTE LA ECUACIÓN DIFERENCIAL

Acudiendo a los datos iniciales sabemos que la producción tiene un máximo de 240 T/año, produciéndose éste a los 25 años.Reflejamos esto matemáticamente como:

dP/dt=0 --------------------------- rLn(k)-rQ'Ln(Q)-rQ'=0 con Q'(25)=240 ya que P era la derivada respecto al tiempo de Q(t)

P_max=240---------------rQ(25)Ln(k/Q(25))=240

Lo que nos proporciona un sistema de ecuaciones con dos incógnitas (r, Q(25))

De este sistema cabe resaltar dos cosas:

1. Al producirse el máximo de la función producción a los 25 años y al depender de la función Q(t),ésta deja de ser una función tomando el valor correspondiente a Q(25).

2. Para calcular la dP/dt debemos utilizar la regla de la cadena, es decir dP/dt=(dP/dQ)*(dQ/dt)

Para resolver este sistema existen múltiples procedimientos, utilizar una calculadora gráfica o la función fsolve de Matlab que lo resuelve iterativamente y de manera aproximada etc

Sin embargo nosotros hemos decidido optar por un programa de Matlab que resuelve ecuaciones mediante cálculo simbólico, obteniendo una solución exacta del sistema. Como observación decir que Matlab no es un programa especializado en cálculo simbólico, si el sistema fuese muy complejo habría que acudir a otro tipo de herramientas, sin embargo en este caso simple e introduciendo las soluciones en las ecuaciones vemos que las verifican sin ningún problema, por lo que decidimos quedarnos con la solución exacta.

Programa Matlab:

syms r q;% Se introducen las variables de manera simbíloca  q=Q(25)
k=10875; 
solve('r*q*ln(10875/q)=240','r*240*ln(10875)-r*240*ln(q)-r*240=0','r','q');
r=ans.r;
r %muestra en pantalla el valor de r
%repito las sentencias por cuestiones de matlab para que me muestre el valor de q
syms r q;
k=10875; 
solve('r*q*ln(10875/q)=240','r*240*ln(10875)-r*240*ln(q)-r*240=0','r','q');
q=ans.q;
q

Lo que nos da unos valores

r=16e/725   que es aproximadamente 0.06

Q(25)=1875/e que es aproximadamente 4000

Conocido el valor de r ya tenemos las expresiones completas tanto de la ecuación diferencial como de la función P(Q)

ANÁLISIS DE LA FUNCIÓN P(Q)

Obtenemos la gráfica de la función P(Q) en Matlab de la siguiente manera

%DATOS
r=(16*exp(1))/725;
k=10875;
Q0=0;
h=0.1;% Introduzco un tamaño de paso inventado suficientemente pequeño para obtener una buena aproximación de la curva
Q=Q0:h:k;% Vector equiespaciado con las toneladas extraibles que será la variable independiente
N=(k-Q0)/h;
P=zeros(1,N+1);%Creo un vector prodcción de ceros que se irá rellenando con la solución
for i=1:N+1;% Bucle que a cada valor del vector Q lo mete en la función producción
    P(i)=r*Q(i)*log(k/Q(i));

  
end

plot(Q,P)%Dibuja la función P(Q)
max(Q) %Muestra en pantalla el máximo de la variable independiente que debe coincidir con k
max(P) %Muestra en pantalla el valor máximo de la función que debe coincidir con 240

De este código resulta la siguiente función

Como era de esperar el máximo lo alcanza en 240T/año para el valor de aproximadamente 4000T que obtuvimos de la resolución del sistema.

Vemos como la función no es simétrica, crece más rápidamente hasta alcanzar el máximo con un decrecimiento posterior más lento, por lo tanto se ajusta bastante bien con las previsiones que esperabamos de la explotación. Aún así se puede comprobar cómo ajusta otro modelo nuestra situación, elegimos por tanto el modelo logístico de Verlhust cuya ecuación diferencial tiene la forma:

dQ/dt=rQ(1-(Q/k))

El tratamiento de la ecuación es análogo al del modelo de Gompertz por lo que no vamos a detallarlo rigurosamente, en esencia debemos obtener las mismas ecuaciones dP/dt=0 (1), 240=rQ(1-(Q/k)) (2) y utilizar el mismo código de Matlab cambiando las ecuaciones. La resolución del sistema nos da unos valores: r*=960/k Q*(25)=k/2 Realizamos un programa en Matlab para comparar ambos resultados y ver cual se ajusta más a nuestra situación:

%verlhust, la implementación es análoga a la del otro método
r2=960/10875; %Nueva tasa de crecimiento
k=10875;
Q0=0;
h=0.1;
Q=Q0:h:k;%Vector variable independiente equiespaciado
N=(k-Q0)/h;
P2=zeros(1,N+1);%Vector de ceros del mismo tamaño que el de la var indep que rellenaré con la solución
for i=1:N+1;%Bucle que asigna a cada valor del vector Q la función Producción (De Verlhust)
    P2(i)=r2*Q(i)*(1-(Q(i)/k));
end
% vuelvo a escribir el modelo de Gompertz para sacar ambos en la misma
% gráfica (Omito explicaciones por ser idéntico al anterior
r=16*exp(1)/725;
k=10875;
Q0=0;
h=0.1;
Q=Q0:h:k;
N=(k-Q0)/h;
P=zeros(1,N+1);
for i=1:N+1;
    P(i)=r*Q(i)*log(k/Q(i));
end
hold on
plot(Q,P2)
plot(Q,P,'r')
legend('Verlhust','Gompertz')

Vemos que el comportamiento según el modelo de Verlhust es muy simétrico. Al esperar una mayor producción durante los primeros años de la explotación es una buena razón para desechar este segundo método.

Podemos hacer una abstracción importante aquí y sacar una buena conclusión, las ecuaciones diferenciales se obtienen en base a principios e hipótesis realizadas, según más acertadas sean nuestras hipótesis mejor reflejará la ecuación el comportamiento real.

Hemos decidido por tanto que el modelo de Gompertz refleja mejor la situación que esperamos en la realidad por lo que a partir de ahora trabajaremos solamente con él.

Una vez hemos obtenido la función P(Q), es muy interesante resolver la ecuación diferencial para conocer la función Q(t). Para ello planteamos el siguiente problema de valor inicial:

dQ/dt=rQLn(k/Q)      Ecuación diferencial

Q(t=0)=0 Condición inicial. Si Q(t) es la cantidad de mineral extraído hasta el tiempo t en el tiempo inicial se habrán extraído 0 T

Ademas de este PVI vamos a añadir una condición que nos viene impuesta por el análisis económico de la explotación. Se sabe que cuando la producción desciende hasta 25 T/año, la explotación deja de ser rentable, Por lo que resolveremos solamente hasta el valor de Q(t) que hace que P(Q)=25, dicho valor de Q lo podemos obtener fácilmente mirando la posición de la columna del vector P que hace a este P=25 y entrar en el vector Q para ese valor de columna.

Este resultado nos da 10450T extraidas que consideraremos como final de la vida útil de la explotación.

Cabe destacar que en la ecuación la función Q(t) que toma el valor 0 en el instante inical, por lo que el cociente se hace infinito y Matlab no lo resolvería, así, cuando implementemos el método aproximaremos el valor de 0 por un valor muy próximo a el, 0.00001.

Para resolver la ecuación vamos a utilizar un método muy rápido y sencillo, el modelo de Euler, este método aproxima la función por la tangente para cada intervalo de tiempo, como nuestro intervalo de tiempo será de un mes, puede valer en primera aproximación pero utilizaremos métodos de mayor orden de error para comprobar la veracidad de este método.

%DATOS INICIALES DEL PVI
r=16*exp(1)/725;
t0=0;
k=10875;
tN=200; %cojo 200 años porque es un tiempo suficientemente largo como para agotar todo el mineral como veremos
N=200*12;%meses que tienen 200 años
h=tN/N; %tamaño de paso de un mes
Q0=0.000001 %No pongo 0 porque al aparecer Q en un cociente daría error el programa
t=t0:h:tN; % vector de tiempos equiespaciado 
Q=zeros(1,N+1);%vector de ceros con lo que será mi solución, y que rellenaré en el bucle
Q(1)=Q0;
for i=1:N %Bucle del método de Euler
    Q(i+1)=Q(i)+h*(r*Q(i)*(log(k)-log(Q(i))));
    if Q(i)>10450 break  %Condición de salida del bucle
    end
end
plot(t,Q)

Con lo que obtenemos la gráfica:

Esta es una función muy interesante y que nos revela mucha información.

En primer lugar llama la atención que durante los primeros 20 años la cantidad de mineral extraído es muy poco, lo que no concuerda con las previsiones que esperábamos ni con la función de producción obtenida que nos decía que la producción crecía durante los primeros 25 años. Por otro lado la función si que se corta(línea vertical) para el valor de 10450 T correspondientes al fin de la vida útil del material

Esto nos hace sospechar de la aproximación de Euler, vamos a comprobarla utilizando métodos mejores como son Heun y Runge Kutta 4 para ver si el problema es del método :

clear all
r=16*exp(1)/725; %Tasa de crecimiento intrínseca
t0=0; %Tiempo inicial
k=10875; %Cantidad máxima extraíble
tN=200; %Tiempo final del intervalo de estudio
N=200*12; %Número de intervalos
h=tN/N; %Anchura de paso 
Q0=0.000001; % Valor inicial
t=t0:h:tN; %Vector de componentes del tiempo
Q=zeros(1,N+1); %Vector de ceros a rellenar por el bucle RG4
Q(1)=Q0; %Introducimos el valor inicial RG4
Q1=zeros(1,N+1); %Vector de ceros a rellenar por el bucle Heun
Q1(1)=Q0; %Introducimos el valor inicial Heun
for i=1:N %Aplicamos el metodo de Runge-Kutta
    k1=r*Q(i)*log(k/Q(i));
    k2=r*(Q(i)+0.5*k1*h)*log(k/(Q(i)+0.5*k1*h));
    k3=r*(Q(i)+0.5*k2*h)*log(k/(Q(i)+0.5*k2*h));
    k4=r*(Q(i)+k3*h)*log(k/(Q(i)+k3*h));
    Q(i+1)=Q(i)+h/6*(k1+2*k2+2*k3+k4);
    if Q(i)>10450, break %Interrumpimos el bucle como consecuencia del final de la vida útil
    end
end
for i=1:N %Aplicamos el metodo de Heun
    k1=r*Q(i)*log(k/Q(i));
    k2=r*(Q(i)+k1*h)*log(k/(Q(i)+k1*h));
    Q1(i+1)=Q(i)+h/2*(k1+k2);
    if Q1(i)>10450 ,break
    end
end
subplot(1,2,1)
plot(t,Q)
title('Grafica Runge-Kutta')
xlabel('Tiempo en años')
ylabel('Toneladas de mineral extraidas')
subplot(1,2,2)
plot(t,Q1)
title('Grafica Heun')
xlabel('Tiempo en años')
ylabel('Toneladas de mineral extraidas')

Podemos observar que las aproximaciones de estos métodos son muy parecidas a la obtenida mediante Euler, por lo que asusmimos que el problema no es del método sino que está en otro lugar.

Por este motivo vamos a mirar el PVI detenidamente:

dQ/dt=rQIn(k/Q)

Q(t=0)=0.000001

Vemos que en la primera iteración este valor tan próximo a 0 de la condición inicial entra multiplicando y diviendo en el cociente al que se le aplica un logarítmo. Dicho cociente se hará muy grande, pero al aplicarle el logaritmo se reducirá bastante y al multiplicar por un valor muy pequeño el resultado se acerca a 0

Por este motivo a la función le cuesta mucho despegarse del 0, vemos que en el momento en que lo hace la función crece rápidamente, por este argumento concretamos que para resolver esta ecuación diferencial es aconsejable utilizar una condición inicial que no se acerque tanto a cero.

Nota: Veremos esto gráficamente al final del texto, por eso no incluimos aquí un ejemplo ilustrativo de esta situación.

Para estudiar la cantidad de mineral extraída cuando el tiempo "t" tiende a infinito asignamos a "t" un valor muy grande (250 años) y observamos que la función tiende a la cantidad máxima extraíble (k=10875),tal y como puede deducirse del modelo de Gompertz, generándose una asíntota horizontal en dicho valor.

clear all
r=16*exp(1)/725;
t0=0;
k=10875;
tN=250; %Valor de t muy grande simulando la tendencia a infinito
N=200*12;
h=tN/N;
Q0=0.000001;
t=t0:h:tN;
Q=zeros(1,N+1);
Q(1)=Q0;
for i=1:N
    k1=r*Q(i)*log(k/Q(i));
    k2=r*(Q(i)+0.5*k1*h)*log(k/(Q(i)+0.5*k1*h));
    k3=r*(Q(i)+0.5*k2*h)*log(k/(Q(i)+0.5*k2*h));
    k4=r*(Q(i)+k3*h)*log(k/(Q(i)+k3*h));
    Q(i+1)=Q(i)+h/6*(k1+2*k2+2*k3+k4);
   
end
plot(t,Q)
max(Q)
title('Aplicacion de RG4 cuando t tiende a infinito')
xlabel('Tiempo en años')
ylabel('Toneladas de mineral extraidas')

Para representar la función de producción respecto del tiempo P(t), obtenemos en primer lugar los valores de la función Q(Cantidad de mineral extraída) respecto del tiempo y aplicamos la función P para estos valores obtenidos.

clear all
r=16*exp(1)/725;
t0=0;
k=10875;
tN=200;
N=200*12;
h=tN/N;
Q0=0.00001;
t=t0:h:tN;
Q=zeros(1,N+1);
Q(1)=Q0;
for i=1:N %Calculamos Q en funcion del tiempo
    Q(i+1)=Q(i)+h*(r*Q(i)*(log(k)-log(Q(i))));
    if Q(i)>10450 ,break
    end
end
P=zeros(1,length(Q));
for i=1:length(Q) %Calculamos P en funcion de los valores de Q obtenidos en el bucle anterior
    P(i)=r*Q(i)*log(k/Q(i));
end
plot(t,P)
title('Gráfica de producción')
xlabel('Tiempo en años')
ylabel('Toneladas por año')

La cantidad que queda sin extraer al terminar la vida util de la explotacion es la diferencia entre la cantidad maxima extraible (10875) y la cantidad extraida hasta el momento (10450).

%OBTENCIÓN DE P(t) UTILIZANDO EL MÉTODO DE HEUN PARA EL MODELO LOGÍSTICO
clear all
%DATOS INICIALES
r=16*exp(1)/725;
t0=0;
k=10875;
tN=500;
N=200*12;
h=tN/N;
Q0=0.00001;
t=t0:h:tN;%Vector equiespaciado de tiempo
Q=zeros(1,N+1);%Vector del mismo tamaño de t que rellenaré con la solución
Q(1)=Q0;%Condición inicial
for i=1:N %Código del método de Heun 
    k1=r*Q(i)*(1-(Q(i)/k));
    k2=r*(Q(i)+k1*h)*(1-(Q(i)+k1*h)/k);
    Q(i+1)=Q(i)+h/2*(k1+k2);
    if Q(i)>10450 break %Para el valor de Q(t) correspondiente a una producción de 25 T/año (vida útil) salgo del bucle
    end
end
%Introducimos los valores de en la función P(Q) para el modelo logístico
P=zeros(1,length(Q));
for i=1:length(Q)
    P(i)=r*Q(i)*log(k/Q(i));
end
plot(t,P)%Dibuja P(t)

clear all
%CAMBIO DEL MÓDELO DEBIDO A LA FALTA DE VERACIDAD DEL MISMO
r=0.069;% Nueva tasa intrínseca de crecimiento
t0=0;
k=9075;%Nueva cantidad máxima extraible
tN=200;
N=200*12;
h=tN/N;
Q0=2695; %Nuevo valor de la condición inicial
t=t0:h:tN;%Vector equiespaciado de tiempos
Q=zeros(1,N+1);%Vector que rellenaré con la solución de Q(t)
Q(1)=Q0;
for i=1:N%Código para el método de Heun para obtener nuevamente Q(t)
    k1=r*Q(i)*log(k/Q(i));
    k2=r*(Q(i)+k1*h)*log(k/(Q(i)+k1*h));
    Q(i+1)=Q(i)+h/2*(k1+k2);
end

subplot(3,1,1)%Junto las gráficas en una misma figura
plot(t,Q,'r')
%P(Q)
Q1=2695:h:k;%Vector que empieza desnde el valor conocido de mineral extraido hasta el agotamiento del mineral
P1=zeros(1,length(Q1));%Vector que rellenaré con la nueva función P(Q)
for i=1:length(Q1) %Bucle que rellena P en con los valores del vector Q
    P1(i)=r*Q1(i)*log(k/Q1(i));
end

%P(Q(t))
P2=zeros(1,length(Q));%Vector que utilizaré para dibujar P(t)
for i=1:length(Q)
    P2(i)=r*Q(i)*log(k/Q(i));
end
%Cómandos para la generación de las figuras en un mismo cuadro para
%facilitar la comparación entre ellas
subplot(3,1,2)
plot(t,P2,'g')
subplot(3,1,3)
plot(Q1,P1,'b')
max(P2)