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

Revisión del 23:11 5 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.

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.

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 si no 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.

Acudiendo a los datos iniciales sabemos que la producción tiene un máximo en 240 T/año, produciendose este a los 25 años.Reflejamos esto matemáticamente coomo:

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

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 e introduciendo las soluciones en las ecuaciones 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)

Vamos a ver la gráfica de la función P(Q) y a analizar resultados

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.

Si análogamente a lo que hemos hecho hasta ahora utilizasemos otra ecuación diferencial para modelizar el problema, como podría ser el modelo de Verlhust cuya ecuación diferencial tiene la forma:

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

obtendríamos otra posible modelización del comportamiento de la explotación. Comparando los dos resultados en una gráfica podemos obtener algunas conclusiones:

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 reflejara 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 con 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 utilizamos el método de Euler con un tamaño de paso equivalente a un mes

%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)

Para dar una mayor precisión a lo obtenido anteriormente se aplican los métodos de Runge-Kutta y Heun. Se puede apreciar que las dos gráficas son completamente idénticas.

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')

Cuando la producción baja a 25 toneladas/año la cantidad de mineral extraído es 10450 toneladas. En este momento la vida útil finaliza, por lo que interrumpimos el bucle con el comando "break".

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)