Diferencia entre revisiones de «Explotación Minera (Grupo 4B)»

Revisión del 15:27 27 abr 2017

1 Introducción

El objetivo de este trabajo es analizar la curva de producción de una explotación minera. Para ello se irán respondiendo a las cuestiones planteadas, analizando los resultados obtenidos y extrayendo conclusiones de dichos resultados.

1.1 ¿Qué ecuación relaciona la producción P y la función Q?

Para poder hallar la relación entre P y Q debemos primero definirlas:

-Q → Cantidad de mineral extraído desde el inicio hasta un tiempo t en toneladas.

-P → Producción en toneladas/año.

Así, la relación entre ambas será la siguiente:

[math] P = \frac{dQ}{dt} [/math]

1.2 Calibrado mediante el modelo de Gompertz y determinación de r

Para poder hallar el coeficiente r debemos tener en cuenta los siguiente datos:

-Ecuación diferencial del modelo de Gompertz

[math] \frac{dQ}{dt} = rQlog\frac{K}{Q} [/math]

Si tenemos en cuenta que hemos definido P como [math] \frac{dQ}{dt} [/math] podemos decir que:

[math] P = rQlog\frac{K}{Q} (1) [/math]

-K es igual a la cantidad total extraíble de mineral, que son 30.800 toneladas.

-La producción máxima es de 510 t/año, por lo que derivando P respecto de Q se haya Q igualando a 0:

[math] \frac{dP}{dQ} = rlog\frac{K}{Q}-r = 0 → r(log\frac{K}{Q}-1) = 0 [/math]

[math] \frac{K}{Q} = e → Q = \frac{K}{e} = 11.331 [/math]

Ahora que tenemos los datos se puede sustituir en la ecuación (1)

[math] 510 = 11.331r → r = 0,045 [/math]

1.3 Obtención de la función P(Q)

La función ya la hemos calculado anteriormente:

[math] P = rQlog\frac{K}{Q} = 0,045Qlog\frac{30.800}{Q} [/math]

Ayudándonos con matlab representamos los valores de P respecto a Q:

%Gráfica Gompertz
t0=0; %tiempo inicial
tN=30800; %tiempo final
q=t0:1:tN;
p1=0.045*q.*log(30800./q); %función definida
plot(q,p1)

centro

Conclusiones de los resultados obtenidos:

Como era de esperar, la producción crece increíblemente rápido cuando la mina se empieza a explotar, sin embargo, al acercarnos a las 11.300 toneladas vemos como la producción alcanza su máximo y desciende drásticamente. Esto se debe a que los recursos están limitados, puesto que la cantidad total extraíble de mineral es de 30.800 toneladas, de tal forma que conforme estos recursos van disminuyendo también lo hace la producción. Al principio la explotación será muy eficiente pero a medida que se consumen los recursos la producción será menor hasta alcanzar un valor de 0 cuando se ha extraído todo el mineral posible (30.800 toneladas).

1.4 Calibrado mediante el modelo de Verhulst

-Ecuación diferencial del modelo de Verhulst

[math] Q' = rQ(1-\frac{Q}{K}) [/math]

Definiendo P como [math] Q' = \frac{dQ}{dt} [/math] podremos decir que:

[math] P = rQ(1-\frac{Q}{K}) (2) [/math]

-K es igual a la cantidad total extraíble de mineral que son 30.800 toneladas.

-La producción máxima es de 510 t/año, y derivando P respecto de Q e igualando a 0 se halla Q:

[math] \frac{dP}{dQ} = r-\frac{2rQ}{K} = 0 → Q = 0,5K = 15.400 t [/math]

Teniendo todos los datos se puede sustituir en la ecuación (2)

[math] 510 = 15.400r(1-\frac{15.400}{30.800}) → r = 0,066 [/math]

Obteniendo como resultado la función P(Q)

[math] P = 0,066Q(1-\frac{Q}{30.800}) [/math]

Trabajando con matlab compararemos las gráficas de ambos resultados:

%Comparación gráfica Gompertz y Verhulst
t0=0; %tiempo inicial
tN=30800; %tiempo final
q=t0:1:tN;
p1=0.045*q.*log(30800./q); %función definida 1
hold on
plot(q,p1)
p2=0.066.*q.*(1-q./30800); %función definida 2
plot(q,p2)
hold off

centro

Conclusiones de los resultados obtenidos:

El descenso de producción según el modelo Verhulst es simétrico respecto al crecimiento, lo que le convierte en un modelo que se adapta peor que el Gompertz, ya que según la demanda de dicho mineral y las características técnicas de su extraccíon se espera un crecimiento rápido de la producción y un descenso más progresivo y lento, condicionado entre otros factores, por debilitamiento de la demanda (también progresivo) y a las dificultades técnicas, cada vez mayores de la extracción.

2 Problema de contorno

La producción comienza a descender a los 27 años (324 meses) del comienzo de la explotación. En ese momento la cantidad de material explotado es Q = 11.332 t, dónde se encuentra al máximo de producción. Luego:

[math] Q(27) = 11.332 t [/math]

Dando lugar al problema de contorno:

[math] PVI=\begin{cases} \frac{dQ}{dt} = 0,045Qlog\frac{30.800}{Q} & \text{}& \\Q(27) = 11.332 & \text{}& \end{cases} [/math]

El tiempo final de la explotación será aquel dónde la producción descienda a 100 t/año, esto es, hasta que se exploten 28490 t. Primero hallamos la función para un tiempo final desconocido:

%Grafico desde t=27 hasra t=?
t0=27;
tN=350;
q0=11332;
f=@(t,q) 0.045*q.*log(30800./q);
h=1/12;
N=round((tN-t0)/h); 
t=t0:h:tN; 
%Definimos el valor q0 inicial 
    q=zeros(size(t));
    q(1)=q0 
    %Euler
 for i=1:N
  q(i+1)=q(i)+h*f(t(i),q(i));
 end
 plot(t,q);
%Grafico desde t=0 hasta t=27
%Lo diferenciamos con el indice 1
t1n=27;
t10=0;
q1n=11332;
t1=t10:h:t1n;
q1=zeros(size(t1));
N1=round((t1n-t10)/h);
q1(N1+1)=q1n;
for i=N1+1:-1:2
    q1(i-1)=q1(i)-h*f(t1(i),q1(i));
end
hold on
plot(t1,q1);
hold off

centro

Ahora debemos saber en qué momento se llega a explotar los 28.490 t, que es a los 84 años.

t0=27;
tN=84;
q0=11332;
f=@(t,q) 0.045*q.*log(30800./q);
h=1/12;
N=round((tN-t0)/h); 
t=t0:h:tN; 
%Definimos el valor q0 inicial 
    q=zeros(size(t));
    q(1)=q0;
    %Euler
 for i=1:N
  q(i+1)=q(i)+h*f(t(i),q(i));
 end
 plot(t,q);
%Grafico desde t=0 hasta t=27
%Lo diferenciamos con el indice 1
t1n=27;
t10=0;
q1n=11332;
t1=t10:h:t1n;
q1=zeros(size(t1));
N1=round((t1n-t10)/h);
q1(N1+1)=q1n;
for i=N1+1:-1:2
    q1(i-1)=q1(i)-h*f(t1(i),q1(i));
end
hold on
plot(t1,q1);
hold off
t(q==28490)

centro

Conclusiones de los resultados obtenidos:

Como se puede apreciar en ambas gráficas los resultados obtenidos son coherentes con el problema planteado, ya que, durante los primeros 27 años la pendiente de la gráfica, que se corresponde con la producción, es muy acusada, tal y cómo se plantea en el problema. La cantidad de material extraído durante los primeros años es muy grande y llegado un momento se estabiliza debido a que los recursos de la cantera son limitados.

2.1 Métodos de Runge Kutta de cuarto orden y Heun

Con las condiciones del apartado anterior nos dedicaremos simplemente a resolver el problema de valor inicial por dos métodos diferentes.

%Runge Kutta
t0=27;
tN=84;
q0=11332;
f=@(t,q) 0.045*q.*log(30800./q);
h=1/12;
N=round((tN-t0)/h); 
t=t0:h:tN; 
%Definimos el valor q0 inicial 
    q=zeros(size(t));
    q(1)=q0
    
    %Runge-Kutta 
    for i=1:N
            K1=f(t(i),q(i));
            K2=f(t(i)+1/2*h,q(i)+1/2*K1*h);
            K3=f(t(i)+1/2*h,q(i)+1/2*K2*h);
            K4=f(t(i)+h,q(i)+K3*h);
            q(i+1)=q(i)+h/6*(K1+2*K2+2*K3+K4);
    end
          hold on
         plot(t,q)
t1n=27;
t10=0;
q1n=11332;
t1=t10:h:t1n;
q1=zeros(size(t1));
N1=round((t1n-t10)/h);
q1(N1+1)=q1n;
for i=N1+1:-1:2
            K1=f(t1(i),q1(i));
            K2=f(t1(i)-1/2*h,q1(i)-1/2*K1*h);
            K3=f(t1(i)-1/2*h,q1(i)-1/2*K2*h);
            K4=f(t1(i)-h,q1(i)-K3*h);
            q1(i-1)=q1(i)-h/6*(K1+2*K2+2*K3+K4);
end
         plot(t1,q1,'b')
         title('Runge Kutta');
         hold off

%Heun
t0=27;
tN=84;
q0=11332;
f=@(t,q) 0.045*q.*log(30800./q);
h=1/12;
N=round((tN-t0)/h); 
t=t0:h:tN; 
%Definimos el valor q0 inicial 
    q=zeros(size(t));
    q(1)=q0
    
    %Heun
          for i=1:N
            K1=f(t(i),q(i));
            K2=f(t(i)+h,q(i)+K1*h);
            q(i+1)=q(i)+h/2*(K1+K2);
          end
          hold on
         plot(t,q)
t1n=27;
t10=0;
q1n=11332;
t1=t10:h:t1n;
q1=zeros(size(t1));
N1=round((t1n-t10)/h);
q1(N1+1)=q1n;
for i=N1+1:-1:2
            K1=f(t1(i),q1(i));
            K2=f(t1(i)-h,q1(i)-K1*h);
            q1(i-1)=q1(i)-h/2*(K1+K2);
          end
         plot(t1,q1,'b')
         title('Heun')
         hold off

Los códigos anteriores nos dan como resultado las gráficas siguientes. Debido al gran parecido entre ambas se ha descartado superponerlas en una misma gráfica.

centro

centro

Conclusiones de los resultados obtenidos:

Ambas gráficas son también a su vez muy similares a la obtenida anteriormente por el método de Euler, no obstante se puede apreciar que en estas dos últimas el crecimiento inicial es mayor, y por lo tanto se ajusta mejor al enunciado del problema. Esto se debe a que tanto el método de Runge Kutta de cuarto orden como el de Heun son de mayor precisión que el de Euler.

2.2 ¿Qué valor límite puede tener Q?

Tenemos que calcular la cantidad de material de la cantera que será extraído pasado un tiempo muy grande. Para ello debemos recordar el problema de contorno anteriormente hallado:

[math] PVI=\begin{cases} \frac{dQ}{dt} = 0,045Qlog\frac{30.800}{Q} & \text{}& \\Q(27) = 11.332 & \text{}& \end{cases} [/math]

Debemos calcular la solución exacta de Q(t) analíticamente. No es necesario ajustar la constante C porque a la hora de calcular el límite es irrelevante. Tras varios cálculos se llega a:

[math] Q(t) = 30.800*e^{-e^{-C-0,045t}} [/math]

Para obtener la cantidad de material obtenido tras un tiempo muy largo debemos calcular el límite de Q(t) cuando t tiende a infinito:

[math] Q_M = \displaystyle\lim_{t \to{+}\infty}{Q(t)} = 30.800*e^{0} = 30.800 t [/math]

En cualquier caso, no habría sido necesario calcularlo analíticamente pues en la gráfica anterior se ve perfectamente que para un tiempo muy grande la cantidad de material extraído tiende a 30.800 toneladas.

Conclusiones de los resultados obtenidos

Se tiene, pues, que la cantidad máxima de material extraído es 30.800 toneladas. El resultado es lógico pues la cantidad total extraíble es de 30.800 toneladas y una vez extraída dicha cantidad los recursos se habrán agotado.

3 Función de producción P(t)

A continuación, analizamos la función de producción P(t), dibujando su curva representativa y obteniendo el punto de máxima producción y los puntos de mayor crecimiento y descenso.Para ello nos basaremos en los datos obtenidos en los apartados anteriores, mediante el método de Euler.

Q=[q1,q];
T=[t1,t];
P= 0.045*Q.*log(30800./Q);
plot(T,P);

centro

Para hallar los puntos de máximo crecimiento y decrecimiento derivamos la función anterior y vemos sus máximos y mínimos, que se corresponderás con los puntos de máximo crecimiento y decrecimiento.

Q=[q1,q];
T=[t1,t];
P2=0.045^2*Q.*log(30800./Q).^2-0.045^2*Q.*log(30800./Q);
plot(T,P2);
T(P2==max(P2))
T(P2==min(P2))

centro

Conclusiones de los resultados obtenidos:

Como era de esperar,la producción crece considerablemente hasta alcanzar un tiempo de 27 años,que es cuando alcanza su máximo, y a partir de ahí la producción comienza a descender. El máximo se alcanza a los 27 años con 509,88 toneladas que es muy próximo a 510(dato proporcionado por el enunciado del problema). El punto de máximo crecimiento se da a los 5 años y 8 meses con una producción de 265,2533 toneladas/año, y el de mínimo se da tras 48 años y 4 meses con una producción de 371,7002 toneladas/año

3.1 Cantidad del mineral sin extraer

Sabemos que la vida útil finaliza cuando la producción es menor de 100 t/año. Previamente hemos hallado la cantidad de mineral extraída hasta ese momento, siendo 28490 t. Si a la cantidad total de mineral extraíble le restamos esta cantidad nos da un valor de 2310 t que quedan sin extraer.

Según el modelo de Verlust la producción alcanza menos de 100 toneladas año cuando se han extraido 29203 toneladas

3.2 Método de Heun

A continuación realizamos la aproximación del modelo anterior mediante el método de Heun.

centro

centro

La gráfica obtenida apenas se diferencia de la anterior perteneciente al método de Euler, no obstante los puntos obtenidos son ligeramente diferentes: el punto de máximo crecimiento se da tras 5 años y 7 meses mientras que el de mínimo tras 48 años y 5 meses.

3.3 Modelo de Verhulst por Heun

t0=27;
tN=84;
q0=11332;
f=@(t,q) 0.066*q.*(1-q/30800);
h=1/12;
N=round((tN-t0)/h); 
t=t0:h:tN; 
%Definimos el valor q0 inicial 
    q=zeros(size(t));
    q(1)=q0;
    
    %Heun
          for i=1:N
            K1=f(t(i),q(i));
            K2=f(t(i)+h,q(i)+K1*h);
            q(i+1)=q(i)+h/2*(K1+K2);
          end
          hold on
         plot(t,q)
t1n=27;
t10=0;
q1n=11332;
t1=t10:h:t1n;
q1=zeros(size(t1));
N1=round((t1n-t10)/h);
q1(N1+1)=q1n;
for i=N1+1:-1:2
            K1=f(t1(i),q1(i));
            K2=f(t1(i)-h,q1(i)-K1*h);
            q1(i-1)=q1(i)-h/2*(K1+K2);
end
         plot(t1,q1,'b')
         title('Heun')
         hold off

Q=[q1,q];
T=[t1,t];
P= 0.066*Q.*(1-Q/30800);
plot(T,P);
max(P)
T(P==max(P))

centro

Q=[q1,q];
T=[t1,t];
P2=0.066^2*Q.*(1-Q/30800)-2*0.066^2*Q.^2/30800.*(1-Q/30800);
plot(T,P2);
T(P2==max(P2))
T(P2==min(P2))

centro

Conclusiones de los resultados obtenidos:

Al aplicar el método de Heun al modelo de Verlhust podemos apreciar como los resultados se alejan de la realidad y el gráfico se desplaza hacia la derecha. Obteniendo el máximo de la producción (508,2 t/año) a los 35 años y 2 meses, cuando debería obtenerse a los 27 años. Los puntos de máximo crecimiento y decrecimiento también se desplazan, obteniéndose respectivamente a los 15 años y 3 meses y a los 55 años y 2 meses. Por tanto, podemos concluir que este modelo no representa fielmente la realidad por lo que no es adecuado. En cuanto a la cantidad de mineral que queda sin extraer, la explotación se parará a los 79 años y 3 meses, habiéndose explotado una cantidad de 29250 toneladas, y quedando por tanto sin explotar una cantidad de 1550 toneladas.

4 Caso final

Se nos propone un último caso en el cuál transcurridos 12 años desde el inicio de la explotación se hace una revisión del comportamiento del modelo, y se ajustan los datos de forma más realista.

Ahora, la cantidad extraída hasta ese año es Q(12) = 8.350 toneladas, y se estima que quedan por extraer en la mina 25.910 toneladas más. Con estos nuevos datos, calculamos los parámetros r y K, y aproximamos de nuevo la función P(t).Obtención de los parámetros r y K.

Para obtener K, sumamos la cantidad de mineral extraída hasta el momento y la cantidad que queda por extraer, y nos queda:

[math] K = 8.350 + 25.910 = 34.260 toneladas [/math]

Para obtener r, se hace un bucle que nos vaya asignando valores de r hasta que la diferencia entre el Q a los 12 años estimado y el real sea muy pequeña.

%valor de r
t0=1;
 tN=12;
 q0=1059.2;%valor inicial obtenido analiticamente de la ley  % de Q(t) tras integrar el PVI y hacer t=0;
 A=3.37;%cte que aparece al integrar Q(t)
 t=t0:1:tN;
 Q=zeros(1,tN);
 Q(1,1)=q0;
 K=34260;
 r=0.01:0.001:1;
 m=length(r);
 for i=1:m
        Q(1,12)=K*exp(-A*exp(-r(i)*t(12)));
  P=8350-Q(1,12);
        if P<=0.01
         break
       end
end
r(i)

De esta forma, tenemos que el valor de r que mejor se adecua a nuestras necesidades ,que es r = 0,073.

A continuación, aproximamos la producción con el nuevo método (método de Heun):

%producción
t0=0;
 tN=100;
 q0=1059.2;
 h=0.01;
 N=(tN-t0)/h;
 t=t0:h:tN;
 q=zeros(1,N+1);
 k=34260;
 r=0.073;
 q(1)=q0;
 p=zeros(1,N+1);
 p(1)=r*q0*log(k/q0);
 for j=1:N
 q(j+1)=q(j)+h*r*q(j)*log(k/q(j));
 end
 for i=1:N
 k1=r^2*q(i)*log(k/q(i))*(log(k/q(i))-1);
 k2=r^2*(q(i)+k1*h)*log(k/(q(i)+k1*h))*(log(k/(q(i)+k1*h))-1);
 p(i+1)=p(i)+h/2*(k1+k2);
 end
 [t',p'];
 figure
 plot(t,p)
 figure
 plot(t,q)
 %Obtenemos el año de máxima producción
 max(p)
 for g=1:N
 if max(p)-p(g)<=0.01
 break
 end
 end
 x=g*tN/N

centro

centro

El punto de máxima producción que obtenemos es P(t) = 920,6 toneladas para el tiempo t = 17 años. Procediendo de forma análoga a la del apartado 9, pero en este caso con los nuevos datos, calculamos la cantidad de mineral extraída hasta el final de la vida útil ,siendo 32.863 toneladas extraídas, quedarán sin extraer por tanto 34.260-32.863 toneladas.

centro

Material sin extraer = 1.397 toneladas

Dibujamos la gráfica simultáneamente de Gompertz con los datos antiguos y nuevos y muestra que con los nuevos datos que nos proporcionan, el modelo tiene un crecimiento y decrecimiento de su producción mucho más rápido que con los valores antiguos, además de que el máximo de producción es mayor y se alcanza antes.

%comparación
clear all
t0=0;
tN=100;%Suponiendo un recorrido representativo de la mina de 100 años
q0=1059.2;%valor  inicial obtenido analiticamente de la ley de  Q(t) tras integrar el PVI
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
q=zeros(1,N+1);
k=30800;
r=0.045;
q(1)=q0;
p=zeros(1,N+1);
p(1)=r*q0*log(k/q0);
for j=1:N %POR  EULER
q(j+1)=q(j)+h*r*q(j)*log(k/q(j));
end
for i=1:N  %POR HEUN
k1=r^2*q(i)*log(k/q(i))*(log(k/q(i))-1);
k2=r^2*(q(i)+k1*h)*log(k/(q(i)+k1*h))*(log(k/(q(i)+k1*h))-1);
p(i+1)=p(i)+h/2*(k1+k2);
end
%valores corregidos
ka=34260;
ra=0.073;
qa(1)=q0;
pa=zeros(1,N+1);
pa(1)=ra*q0*log(ka/q0);
for j=1:N
qa(j+1)=qa(j)+h*ra*qa(j)*log(ka/qa(j));
end
for i=1:N
k1=ra^2*qa(i)*log(ka/qa(i))*(log(ka/qa(i))-1);
k2=ra^2*(qa(i)+k1*h)*log(ka/(qa(i)+k1*h))*(log(ka/(qa(i)+k1*h))-1);
pa(i+1)=pa(i)+h/2*(k1+k2);
end
hold on 
plot(t,p)
plot(t,pa,'r')
legend ('modelo antiguo','modelo corregido')
hold off

centro