Diferencia entre revisiones de «Explotación minera (Grupo 5C)»

Revisión del 12:25 5 mar 2015

1 Introducción

• El problema nos pide el análisis de la explotación de un yacimiento de mineral. Dicha explotación sigue un modelo logístico de Gompertz, cuya ecuación tiene la siguiente forma:
  

:[math] {dQ \over dt}=rQlog ({K \over Q}) [/math]

donde Q(t) es la cantidad de mineral extraído, K la cantidad total extraíble y r la tasa intrínseca de crecimiento.

• En nuestro caso, sabemos que la extracción de mineral tendrá un crecimiento muy rápido de producción durante los primeros 25 años, momento a partir del cual descenderá lentamente debido a diversos factores. Además de esto, conocemos la cantidad total extraíble del yacimiento, por lo que nuestra ecuación inicial quedará de la siguiente forma:
  

:[math] {dQ \over dt}=rQlog ({10875 \over Q}) [/math]

2 Relación entre cantidad y producción

Mediante las siguientes gráficas, se demuestra que la relación entre Q y P es la siguiente:

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

o lo que es lo mismo:[math] P=Q' [/math]

• 

  Descripción1

• 

  Descripción2

3 Modelo de Gompertz: Valor de la tasa de crecimiento (r) y resolución de la función de producción

3.1 Tasa de crecimiento r

• Partimos del siguiente dato inicial:

K= 10875 toneladas

• La ecuación del modelo de Gompertz es:
  

:[math]{dQ \over dt}=rQlog({K \over Q})[/math]

• Condición:

P(t=25 años)= 240 t/año [math] \rightarrow [/math] P(25)=240 [math] \rightarrow [/math] Q'(25)=240

• Resolución:

La solución de la ecuación de Gompertz es la siguiente:

: [math] Q(t)=Ke^{-Ae^{-rt}} [/math]

El enunciado dice que la producción es máxima a los 25 años, por lo que si Q' tiene un máximo en un punto, Q"=0 en dicho punto.
Derivamos la ecuación inicial del modelo y la igualamos con la condición que hemos obtenido anteriormente.

: [math] Q"=r^2Qln({K \over Q})ln({K \over Q}-1)=0  [/math]

Despejando esta ecuación llegamos a que Q(25 años)= 4000,69 t.
Conociendo Q(25) podemos despejar "r" de nuestro problema:

:  [math] Q'= rQ(25)ln({10375 \over Q(25)})=240 [/math]

Así obtenemos r=0,05999 [math] \simeq [/math] 0,06
Con estos datos, ya podemos despejar A de [math] Q(t)=Ke^{-Ae^{-rt}} [/math], siendo A=4,4705
Por lo tanto, nuestra ecuación del modelo de Gompertz quedará de la siguiente forma:

:  [math] Q(t)=10875e^{-4,47e^{-0,6t}} [/math]

3.2 Función de producción P(Q)

Para obtener la función P(Q), obtendremos por un lado P(t) y por otro Q(t), para después representarlos en la forma P(Q).

• 

  Descripción1

CODIGO MATLAB

%Lo que haremos será calcular valores de Q para distintos valores de t, y
%luego calcular valores de P para esos mismo valores de t. Una vez tengamos
%los valores de P y Q en sendos vectores, los representaremos
K=10875;
A=4.47;
r=0.06;
t=0:1:100;
Q=zeros(1,101);
P=zeros(1,101);
for i=1:101
Q(1,i)=K*exp(-A*exp(-r*t(i)));
end
for j=1:101
P(j)=r*Q(j)*log(K/Q(j));
end
plot(Q,P,'*g')
xlabel('Cantidad de mineral extraída en toneladas')
ylabel('Producción de mineral en toneladas/año')

4 Modelo de Verhlust: Valor de la tasa de crecimiento (r) y resolución de la función de producción

4.1 Tasa de crecimiento r

• Partimos del siguiente dato inicial:

K= 10875 toneladas

• La ecuación del modelo de Verhlust es:
  

[math]{dQ \over dt}=rQ(1-{Q \over K})[/math]

• Condición:

P(t=25 años)= 240 t/año [math] \rightarrow [/math] P(25)=240 [math] \rightarrow [/math] Q'(25)=240

• Resolución:

Para la resolución, utilizamos el mismo planteamiento que en el apartado anterior del modelo de Gompertz.

: [math] Q"=rQ'(1-{Q \over 5437,5})=0  [/math]

En el caso del modelo de Verhlust obtenemos r=0,0883 y Q(25 años)=5437,5 t. Así, la ecuación del modelo tendrá la siguiente forma:

: [math] Q'=0,0883Q(1-{Q \over 10875})  [/math]

4.2 Función de producción P(Q)

Como sabemos que P=Q', la ecuación de P(t) es la siguiente:

: [math] P(t)=0,0883Q(t)(1-{Q(t) \over 10875})=0  [/math]

La ecuación de Q(t) la hemos obtenido integrando la del modelo obtenida anteriormente, y tiene la forma:

: [math] Q(t)={1196,25e^0,088276t \over 1+11e^0.088276t}  [/math]

Al igual que en el apartado anterior, obtenemos P(Q) dibujando ambas funciones en la misma gráfica:

• 

  Descripción1

[CODIGO MATLAB]

5 Comparación de P(Q) según modelos de Gompertz y Verhlust

6 Problema de valor inicial utilizando distintos métodos

6.1 Método de Euler

CODIGO MATLAB

%Apartado 5. Metodo de Euler
r=0.06;
k=10875;
t0=0;
tN=100; %Consideramos periodo de representacion de 100 años como suficientemente representativo graficamente
h=(1/12); %Nos piden que h sea de 1 mes. Hay 1/12 meses en un año
t=t0:h:tN; %t tiene 1201 elementos. 1201 eltos - 1 elemento inicial = 1200 eltos. 1200 elementos entre 100 años que hemos representado = 12 eltos.
%12 elementos por año, es decir, uno por cada mes. OK!!
Q0=124.4896; %Suponemos para Euler el mismo valor inicial que el obtenido con la condicion a t=25 años
%del problema sustituida en la ec. Gompertz programada segun la sol. analitica
Q=zeros(1,length(t));
P=zeros(1,length(t));
Q(1,1)=Q0;
P(1,1)=r*Q(1)*log(k/Q(1));
for i=1:length(t)-1
    Q(1,i+1)=Q(1,i)+h*r*Q(1,i)*log(k/Q(1,i));
    P(1,i+1)=r*Q(1,i+1)*log(k/Q(1,i+1));
end
for n=1:length(t)-1
    P(1,n);
    if P(1,n)>=25
        Plim=P(1,n);
        m=n;
    end
end
Plim %Plim es el ultimo valor de la produccion que esta por encima de las 25t/año
m %Paso correspondiente al valor de Plim, se puede comprobar que P(m)=Plim=25t/año

%Ahora que ya conocemos el paso de tiempo en el que finaliza la vida util
%de la explotacion (m=787), se puede reescribir la funcion Q(t), que valdra
%0 a partir de dicho paso de tiempo
for j=m:length(t)-1
    Q(1,j+1)=0;
end

figure
plot(t,Q(1,:))
xlabel('Tiempo en años')
ylabel('Cantidad de material extraido en toneladas')
legend('Cantidad de material extraido (Gompertz / Euler)','Location','Best')
figure
hold on
plot(t,P(1,:))
plot(t,t.*0+25,'r')
xlabel('Tiempo en años')
ylabel('Produccion en toneladas/año')
legend('Produccion (Gompertz / Euler)','Produccion limite (por debajo la vida util de la explotacion habra terminado)','Location','Best')
hold off

6.2 Método de Runge-Kutta

%Apartado 6. Metodos de Runge Kutta y Heun aplicados al apartado anterior clear all r=0.06; k=10875; t0=0; tN=100; %Consideramos periodo de representacion de 100 años como suficientemente representativo graficamente h=(1/12); %Nos piden que h sea de 1 mes. Hay 1/12 meses en un año t=t0:h:tN; %t tiene 1201 elementos. 1201 eltos - 1 elemento inicial = 1200 eltos. 1200 elementos entre 100 años que hemos representado = 12 eltos. 12 elementos por año, es decir, uno por cada mes. OK!! Q0=124.4896; %Suponemos para Euler el mismo valor inicial que el obtenido con la condicion a t=25 años del problema sustituida en la ec. Gompertz programada segun la sol. analitica

%Runge Kutta Q_rk=zeros(1,length(t)); P_rk=zeros(1,length(t)); Q_rk(1,1)=Q0; P_rk(1,1)=r*Q_rk(1)*log(k/Q_rk(1));

%Heun Q_heun=zeros(1,length(t)); P_heun=zeros(1,length(t)); Q_heun(1,1)=Q0; P_heun(1,1)=r*Q_heun(1)*log(k/Q_heun(1));

%Aproximacion por Runge Kutta for i=1:length(t)-1

   K1=r*Q_rk(1,i)*log(k/Q_rk(1,i));
   K2=r*(Q_rk(1,i)+0.5*K1*h)*log(k/(Q_rk(1,i)+0.5*K1*h));
   K3=r*(Q_rk(1,i)+0.5*K2*h)*log(k/(Q_rk(1,i)+0.5*K2*h));
   K4=r*(Q_rk(1,i)+K3*h)*log(k/(Q_rk(1,i)+K3*h));
   Q_rk(1,i+1)=Q_rk(1,i)+h/6*(K1+2*K2+2*K3+K4);
   P_rk(1,i+1)=r*Q_rk(1,i+1)*log(k/Q_rk(1,i+1));

end for n=1:length(t)-1

   P_rk(1,n);
   if P_rk(1,n)>=25
       Plim=P_rk(1,n);
       m=n;
   end

end

for j=m:length(t)-1

   Q_rk(1,j+1)=0;

end

%Aproximacion por Heun for i=1:length(t)-1

   k1=r*Q_heun(1,i)*log(k/Q_heun(1,i));
   k2=r*(Q_heun(1,i)+k1*h)*log(k/(Q_heun(1,i)+k1*h));
   Q_heun(1,i+1)=Q_heun(1,i)+h/2*(k1+k2);
   P_heun(1,i+1)=r*Q_heun(1,i+1)*log(k/Q_heun(1,i+1));

end for n=1:length(t)-1

   P_heun(1,n);
   if P_heun(1,n)>=25
       Plim=P_heun(1,n);
       m=n;
   end

end for j=m:length(t)-1

   Q_heun(1,j+1)=0;

end

%Representacion Runge Kutta figure plot(t,Q_rk(1,:),'c') xlabel('Tiempo en años') ylabel('Cantidad de material extraido en toneladas') legend('Cantidad extraida (Gompertz / Runge Kutta)','Location','Best') figure hold on plot(t,P_rk(1,:),'c') plot(t,t.*0+25,'r') xlabel('Tiempo en años') ylabel('Produccion en toneladas/año') legend('Produccion (Gompertz / Runge Kutta)','Produccion limite (por debajo la vida util de la explotacion habra terminado)','Location','Best') hold off

%Representacion Heun figure plot(t,Q_heun(1,:),'k') xlabel('Tiempo en años') ylabel('Cantidad de material extraido en toneladas') legend('Cantidad extraida (Gompertz / Heun)','Location','Best') figure hold on plot(t,P_heun(1,:),'k') plot(t,t.*0+25,'r') xlabel('Tiempo en años') ylabel('Produccion en toneladas/año') legend('Produccion (Gompertz / Heun)','Produccion limite (por debajo la vida util de la explotacion habra terminado)','Location','Best') hold off

%Representacion simultanea de Runge Kutta y Heun figure hold on plot(t,Q_rk(1,:),'c') plot(t,Q_heun(1,:),'k') xlabel('Tiempo en años') ylabel('Cantidad de material extraido en toneladas') legend('Cantidad extraida (Gompertz / Runge Kutta)','Cantidad extraida (Gompertz / Heun)','Location','Best') hold off %Hay que hacer muchisimo zoom para empezar a ver diferencas entre las dos %funciones!!! figure hold on plot(t,P_rk(1,:),'c') plot(t,P_heun(1,:),'k') plot(t,t.*0+25,'r') xlabel('Tiempo en años') ylabel('Produccion en toneladas/año') legend('Produccion (Gompertz / Runge Kutta)','Produccion (Gompertz / Heun)','Produccion limite (por debajo la vida util de la explotacion habra terminado)','Location','Best')

6.3 Método de Heun

7 Función de producción P(t): Modelo de Gompertz

• Utilizando el método de Heun aproximamos la función de producción P(t).
  
• A continuación obtenemos el punto de máxima producción, que es P(t)=240,18 para t=24,18 años.
  
• Para los puntos de mayor crecimiento y descenso hacemos...
  

• 

  Representación de la función P(t)

• 

  Comparación de P(t) obtenida analiticamente y de forma numérica

CODIGO MATLAB

clear all
t0=0;
tN=100;
q0=124.5;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
q=zeros(1,N+1);
k=10875;
r=0.06;
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'];
plot(t,p)
%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
y=g*tN/N;

7.1 Cantidad de mineral que queda por extraer

8 Función de producción P(t): Modelo de Verhlust

• En este caso volvemos a utilizar el método de Heun para aproximar la función.
  
• El punto de máxima producción que obtenemos es P(t)=240,18 para un tiempo t=50,38 años
  
• Los puntos de máximo y crecimiento son ???
  

• 

  Representación de P(t)

• 

  Comparación de P(t) obtenida de forma analítica y aproximada numéricamente

CODIGO MATLAB:

clear all
t0=0;
tN=100;
q0=124.5;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
q=zeros(1,N+1);
k=10875;
r=0.0883;
q(1)=q0;
p=zeros(1,N+1);
p(1)=r*q0*(1-(q0/k));
for j=1:N
q(j+1)=q(j)+h*r*q(j)*(1-(q(j)/k));
end
for i=1:N
k1=r*p(i)*(1-(2*q(i)/k));
k2=r*(p(i)*(1-(2*(q(i)+k1*h)/k)));
p(i+1)=p(i)+h/2*(k1+k2);
end
[t',p'];
plot(t,p)
%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
y=g*tN/N;

8.1 Cantidad de mineral que queda sin extraer

9 Comparación de P(t) según Gompertz y Verhlust

10 Aproximación del modelo con nuevos datos

• El enunciado nos dice que transcurridos 12 años del inicio de la explotación, se hace una revisión del comportamiento del modelo y se ajustan los datos de forma más real.
  
• Ahora, la cantidad extraída hasta ese año 12 es Q(12)=2695 toneladas, y se estima que quedan por extraer en la mina 9075 toneladas más.
  
• Con estos nuevos datos, calculamos los parámetros r y K, y aproximamos la función P(t).

10.1 Obtención de los parámetros r y K

• Para obtener K, simplemente sumamos la cantidad de mineral extraída hasta el momento y la cantidad que queda por extraer, y nos queda [math]K=2695+9075=11770 [/math]toneladas.

10.2 Producción

• Volvemos a aproximar el modelo utilizando el método de Heun.
  
• El punto de máxima producción que obtenemos es P(t)=400,87 toneladas para el tiempo t=16,33 años.
  
• La cantidad de mineral que queda sin extraer la obtenemos...???

• 

  Representación de P(t)

CODIGO MATLAB

clear all
t0=0;
tN=100;
q0=124.5;
h=0.01;
N=(tN-t0)/h;
t=t0:h:tN;
q=zeros(1,N+1);
k=11770;
r=0.0925;
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;

10.3 Comparación del nuevo modelo con la previsión inicial

• 

  Gráfica comparativa