[[ Check out my Wordpress blog Context/Earth for environmental and energy topics tied together in a semantic web framework ]]

Thursday, November 03, 2005

General Shock Model Program

I attached source code for a generalized command-line driven program that takes input discovery data (+ other parameters) and which then generates production curves from the date of first discovery to the present time (and on out into the future). The code basically refactors the specific code from other posts in the Oil Shock Model series, to get a better separation of data from the algorithms.

The program reads the input parameters from either a set of files (for yearly data) or from the term following the option (for fixed constants). It doesn't matter what units you use as long as they remain consistent with the other ones. The program sends the output data to the shell terminal.

File containing discovery data (continuous years, one entry per line)

File containing shock data (only years with discontinuities), two entries per line, "year" + (space) + "extraction rate that year". The program interpolates between entries, so to prevent shocks, give it only two lines (start and finish times) with the same extraction rate on each line.

Optional file containing production data (continuous years, one entry per line). The program doesn't do anything with this data; it only carries it over to make it convenient for comparison to the estimated production.

The 1/e rate for remaining fallow (i.e. 1 over the average time between discovery and construction)

The 1/e rate for finishing construction (i.e. 1 over the build time)

The 1/e rate for reaching maturing (i.e. 1 over the startup time)

Timespan for complete simulation

Typical invocation:
shockmodel.exe -1 0.2 -2 0.2 -3 0.2 -t 600 -d strikes.dat -p prod.dat -s shocks.dat > output

with Text_IO;
with Ada.Numerics.Elementary_Functions;
with GNAT.Command_Line;
with Ada.Float_Text_IO;

procedure ShockModel is
type Flt is array (Natural range <>) of Float;
type Flt_Access is access all Flt;

Expansion_Factor : constant := 10; -- Year subtics

Time_Span : Integer := 600;

-- Safe array retrieval function
function Get (Arr : Flt; I : Integer) return Float is
if I in Arr'Range then
return Arr (I);
return 0.0;
end if;
end Get;

-- Discretized convolution function
function "*" (A, B : in Flt) return Flt is
Total : constant Natural := A'Length + B'Length;
C : Flt (0 .. Total);
V : Float;
for J in 0 .. Total loop
V := 0.0;
for I in 0 .. J loop
V := V + Get (A, I) * Get (B, J - I);
end loop;
C (J) := V;
end loop;
return C;
end "*";

-- Discretized exponential function
function Exponential (L : Natural; Alpha : Float) return Flt is
use Ada.Numerics.Elementary_Functions;
R : Flt (0 .. L - 1);
for I in 0 .. L - 1 loop
R (I) := Alpha * Exp (-Alpha * Float (I));
end loop;
return R;
end Exponential;

-- Fills in data between discrete years
function Expander (Arr : Flt) return Flt is
R : Flt (0 .. Arr'Length * Expansion_Factor - 1);
for I in R'Range loop
R (I) := Arr (I / Expansion_Factor);
end loop;
return R;
end Expander;

-- Safe function for reading forcing function
function Discovery_Window
(Discovery : Flt;
Year : Float;
Start : Float)
return Float
Y : constant Float := Year - Start;
if Y < 0.0 then
return 0.0;
end if;
return Discovery (Integer (Y) * Expansion_Factor); -- Subdivide by factor
end Discovery_Window;

type Shock is record
Year : Float; -- Year of shock
Rate : Float; -- Rate of extraction
end record;
type Shocks is array (Natural range <>) of Shock;
type Shocks_Access is access all Shocks;

-- Make the shock a continuous function
-- by interpolating between points in the list
function Interpolate_Shocks (Year : Float; S : Shocks) return Float is
K : Integer := S'Last - 1;
for J in S'First + 1 .. S'Last loop
if S (J).Year > Year then
K := J - 1;
end if;
end loop;
return S (K).Rate +
(S (K + 1).Rate - S (K).Rate) * (Year - S (K).Year) /
(S (K + 1).Year - S (K).Year);
end Interpolate_Shocks;

-- Strikes into Mature Discovery
Strikes : Flt_Access;
Discovery : Flt_Access;

Phases : Flt (1 .. 3) := (others => 0.0);

-- Generates the mature discovery by repeated convolution
procedure Phasing is
Avg_Fallow : constant Flt :=
Exponential (Time_Span, Phases (1) / Float (Expansion_Factor));
Avg_Build : constant Flt :=
Exponential (Time_Span, Phases (2) / Float (Expansion_Factor));
Avg_Mature : constant Flt :=
Exponential (Time_Span, Phases (3) / Float (Expansion_Factor));
Discovery :=
new Flt'(Expander (Strikes.all) * Avg_Fallow * Avg_Build * Avg_Mature);
end Phasing;

S : Shocks_Access;
Production : Flt_Access; -- To compare against, optional

procedure Compute is
C : Float := 0.0; -- Reserve
Rate : Float; -- Extraction Rate
P : Float; -- Production Rate
DT : constant Float := 0.01; -- 1/100th of year
Start : constant Float := S (S'First).Year;
Finish : constant Float := S (S'Last).Year;
Time : Float := Start;
V : Float := 0.0; -- Volume of oil extracted
Data : Float := 0.0;
PIndex : Integer;
-- Integration of discovery window against extraction rate
Rate := Interpolate_Shocks (Time, S.all);
C := C +
Discovery_Window (Discovery.all, Time, Start) * DT -
C * Rate * DT;
P := Rate * C;
V := V + Rate * C * DT;
if Time > Float (Integer (Time)) - DT / 2.0 and
Time < Float (Integer (Time)) + DT / 2.0
then -- Print output only once per year
PIndex := Integer (Time - Start);
if PIndex in Production'Range then
Data := Production (PIndex);
Data := 0.0;
end if;
(Integer'Image (Integer (Time)) & "," & P'Img & "," & Data'Img);
end if;
Time := Time + DT;
exit when Time > Finish;
end loop;
end Compute;

function Get_Flt_Data (File : in String) return Flt_Access is
FT : Text_IO.File_Type;
function Get_Data (Data : in Flt) return Flt is
Value : Flt (0 .. 0);
Ada.Float_Text_IO.Get (FT, Value (0));
Text_IO.Skip_Line (FT);
return Get_Data (Data & Value);
when Text_IO.End_Error =>
return Data;
end Get_Data;
Text_IO.Open (FT, Text_IO.In_File, File);
return new Flt'(Get_Data (Flt'(1 .. 0 => 0.0)));
end Get_Flt_Data;

function Get_Shock_Data (File : in String) return Shocks_Access is
FT : Text_IO.File_Type;
function Get_Data (Data : in Shocks) return Shocks is
Pair : Shocks (0 .. 0);
Ada.Float_Text_IO.Get (FT, Pair (0).Year);
Ada.Float_Text_IO.Get (FT, Pair (0).Rate);
Text_IO.Skip_Line (FT);
return Get_Data (Data & Pair);
when Text_IO.End_Error =>
return Data;
end Get_Data;
Text_IO.Open (FT, Text_IO.In_File, File);
return new Shocks'(Get_Data (Shocks'(1 .. 0 => (0.0, 0.0))));
end Get_Shock_Data;

-- Command line options
case GNAT.Command_Line.Getopt ("1: 2: 3: t: d: s: p:") is
when ASCII.NUL =>
when '1' =>
Phases (1) := Float'Value (GNAT.Command_Line.Parameter);
when '2' =>
Phases (2) := Float'Value (GNAT.Command_Line.Parameter);
when '3' =>
Phases (3) := Float'Value (GNAT.Command_Line.Parameter);
when 't' =>
Time_Span := Integer'Value (GNAT.Command_Line.Parameter);
when 'd' =>
Strikes := Get_Flt_Data (GNAT.Command_Line.Parameter);
when 's' =>
S := Get_Shock_Data (GNAT.Command_Line.Parameter);
when 'p' =>
Production := Get_Flt_Data (GNAT.Command_Line.Parameter);
when others =>
raise Program_Error; -- cannot occur!
end case;
end loop;
when GNAT.Command_Line.Invalid_Switch =>
Text_IO.Put_Line ("Invalid Switch " & GNAT.Command_Line.Full_Switch);
when GNAT.Command_Line.Invalid_Parameter =>
Text_IO.Put_Line ("No parameter for " & GNAT.Command_Line.Full_Switch);
end ShockModel;


Professor Anonymous Anonymous said...

Sweet! thank you!

4:55 PM  
Professor Blogger @whut said...

You are most welcome.

5:13 PM  
Professor Blogger Khebab said...

I've implemented your code in Matlab and R. I've tried to validate my implementation using your example on ASPO data:


I retrieved your result only if:
Expansion_Factor= 1;
if Expansion_Factor= 10, the shock model is giving smaller values by 10% compared to the real data.

8:27 AM  
Professor Blogger @whut said...

The "expansion factor" is only used to get better time resolution for high rates. It's a kind of Nyquist factor hack for subdividing a data set into subyear increments.

The best way to validate the porting of the code is to take a delta discovery value (K*delta(t)) and make all the rates the same (the 3 stages plus the depletion rate). The resultant curve should lay exactly on top of a Gamma distribution of order N=4. This is essentially an exponential convolved 3 times with itself. This should show an analytical equivalence.

I mention this because it well tell you if something got screwed up in the expansion factor algorithm code translation. Otherwise, keeping the factor at 1 is OK as long as the numeric errors don't creep in as you increase the rates.

BTW, I also have Python code for this algorithm.

3:45 PM  

Post a Comment

<< Home

"Like strange bulldogs sniffing each other's butts, you could sense wariness from both sides"